overlap.F90

!-*- mode: F90 -*-!
!------------------------------------------------------------!
! This file is distributed as part of the Wannier90 code and !
! under the terms of the GNU General Public License. See the !
! file `LICENSE' in the root directory of the Wannier90      !
! distribution, or http://www.gnu.org/copyleft/gpl.txt       !
!                                                            !
! The webpage of the Wannier90 code is www.wannier.org       !
!                                                            !
! The Wannier90 code is hosted on GitHub:                    !
!                                                            !
! https://github.com/wannier-developers/wannier90            !
!------------------------------------------------------------!

module w90_overlap
  !! This module reads in the overlap (Mmn) and Projections (Amn)
  !! and performs simple operations on them.

  use w90_constants, only: dp, cmplx_0, cmplx_1
  use w90_parameters, only: disentanglement
  use w90_io, only: stdout
  use w90_comms, only: on_root, comms_bcast

  implicit none

  private

!~  public :: overlap_dis_read
  public :: overlap_allocate
  public :: overlap_read
  public :: overlap_dealloc
  public :: overlap_project
  public :: overlap_project_gamma  ![ysl]
!~  public :: overlap_check_m_symmetry

contains

  !%%%%%%%%%%%%%%%%%%%%%
  subroutine overlap_allocate()
    !%%%%%%%%%%%%%%%%%%%%%
    !! Allocate memory to read Mmn and Amn from files
    !! This must be called before calling overlap_read

    use w90_parameters, only: num_bands, num_wann, num_kpts, nntot, timing_level, &
      u_matrix, m_matrix_orig, m_matrix_orig_local, a_matrix, &
      u_matrix_opt, m_matrix, m_matrix_local
    use w90_io, only: io_error, io_stopwatch
    use w90_comms, only: my_node_id, num_nodes, comms_array_split

    integer :: ierr
    ! Needed to split an array on different nodes
    integer, dimension(0:num_nodes - 1) :: counts
    integer, dimension(0:num_nodes - 1) :: displs

    if (timing_level > 0) call io_stopwatch('overlap: allocate', 1)

    call comms_array_split(num_kpts, counts, displs)

    allocate (u_matrix(num_wann, num_wann, num_kpts), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating u_matrix in overlap_read')
    u_matrix = cmplx_0

    if (disentanglement) then
      if (on_root) then
        allocate (m_matrix_orig(num_bands, num_bands, nntot, num_kpts), stat=ierr)
        if (ierr /= 0) call io_error('Error in allocating m_matrix_orig in overlap_read')
      endif
      allocate (m_matrix_orig_local(num_bands, num_bands, nntot, counts(my_node_id)), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating m_matrix_orig_local in overlap_read')
      allocate (a_matrix(num_bands, num_wann, num_kpts), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating a_matrix in overlap_read')
      allocate (u_matrix_opt(num_bands, num_wann, num_kpts), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating u_matrix_opt in overlap_read')
    else
      if (on_root) then
        allocate (m_matrix(num_wann, num_wann, nntot, num_kpts), stat=ierr)
        if (ierr /= 0) call io_error('Error in allocating m_matrix in overlap_read')
        m_matrix = cmplx_0
      endif
      allocate (m_matrix_local(num_wann, num_wann, nntot, counts(my_node_id)), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating m_matrix_local in overlap_read')
      m_matrix_local = cmplx_0
    endif

    if (timing_level > 0) call io_stopwatch('overlap: allocate', 2)

  end subroutine overlap_allocate

  !%%%%%%%%%%%%%%%%%%%%%
  subroutine overlap_read()
    !%%%%%%%%%%%%%%%%%%%%%
    !! Read the Mmn and Amn from files
    !! Note: one needs to call overlap_allocate first!

    use w90_parameters, only: num_bands, num_wann, num_kpts, nntot, nncell, nnlist, &
      num_proj, lselproj, proj2wann_map, &
      devel_flag, u_matrix, m_matrix, a_matrix, timing_level, &
      m_matrix_orig, u_matrix_opt, cp_pp, use_bloch_phases, gamma_only, & ![ysl]
      m_matrix_local, m_matrix_orig_local, lhasproj
    use w90_io, only: io_file_unit, io_error, seedname, io_stopwatch
    use w90_comms, only: my_node_id, num_nodes, &
      comms_array_split, comms_scatterv

    implicit none

    integer :: nkp, nkp2, inn, nn, n, m, i, j
    integer :: mmn_in, amn_in, num_mmn, num_amn
    integer :: nnl, nnm, nnn, ncount
    integer :: nb_tmp, nkp_tmp, nntot_tmp, np_tmp, ierr
    real(kind=dp) :: m_real, m_imag, a_real, a_imag, mu_tmp, sigma_tmp
    complex(kind=dp), allocatable :: mmn_tmp(:, :)
    character(len=50) :: dummy
    logical :: nn_found
    ! Needed to split an array on different nodes
    integer, dimension(0:num_nodes - 1) :: counts
    integer, dimension(0:num_nodes - 1) :: displs

    if (timing_level > 0) call io_stopwatch('overlap: read', 1)

    call comms_array_split(num_kpts, counts, displs)

    if (disentanglement) then
      if (on_root) then
        m_matrix_orig = cmplx_0
      endif
      m_matrix_orig_local = cmplx_0
      a_matrix = cmplx_0
      u_matrix_opt = cmplx_0
    endif

    if (on_root) then

      ! Read M_matrix_orig from file
      mmn_in = io_file_unit()
      open (unit=mmn_in, file=trim(seedname)//'.mmn', &
            form='formatted', status='old', action='read', err=101)

      if (on_root) write (stdout, '(/a)', advance='no') ' Reading overlaps from '//trim(seedname)//'.mmn    : '

      ! Read the comment line
      read (mmn_in, '(a)', err=103, end=103) dummy
      if (on_root) write (stdout, '(a)') trim(dummy)

      ! Read the number of bands, k-points and nearest neighbours
      read (mmn_in, *, err=103, end=103) nb_tmp, nkp_tmp, nntot_tmp

      ! Checks
      if (nb_tmp .ne. num_bands) &
        call io_error(trim(seedname)//'.mmn has not the right number of bands')
      if (nkp_tmp .ne. num_kpts) &
        call io_error(trim(seedname)//'.mmn has not the right number of k-points')
      if (nntot_tmp .ne. nntot) &
        call io_error(trim(seedname)//'.mmn has not the right number of nearest neighbours')

      ! Read the overlaps
      num_mmn = num_kpts*nntot
      allocate (mmn_tmp(num_bands, num_bands), stat=ierr)
      if (ierr /= 0) call io_error('Error in allocating mmn_tmp in overlap_read')
      do ncount = 1, num_mmn
        read (mmn_in, *, err=103, end=103) nkp, nkp2, nnl, nnm, nnn
        do n = 1, num_bands
          do m = 1, num_bands
            read (mmn_in, *, err=103, end=103) m_real, m_imag
            mmn_tmp(m, n) = cmplx(m_real, m_imag, kind=dp)
          enddo
        enddo
        nn = 0
        nn_found = .false.
        do inn = 1, nntot
          if ((nkp2 .eq. nnlist(nkp, inn)) .and. &
              (nnl .eq. nncell(1, nkp, inn)) .and. &
              (nnm .eq. nncell(2, nkp, inn)) .and. &
              (nnn .eq. nncell(3, nkp, inn))) then
            if (.not. nn_found) then
              nn_found = .true.
              nn = inn
            else
              call io_error('Error reading '//trim(seedname)// &
                            '.mmn. More than one matching nearest neighbour found')
            endif
          endif
        end do
        if (nn .eq. 0) then
          if (on_root) write (stdout, '(/a,i8,2i5,i4,2x,3i3)') &
            ' Error reading '//trim(seedname)//'.mmn:', ncount, nkp, nkp2, nn, nnl, nnm, nnn
          call io_error('Neighbour not found')
        end if
        if (disentanglement) then
          m_matrix_orig(:, :, nn, nkp) = mmn_tmp(:, :)
        else
          ! disentanglement=.false. means numbands=numwann, so no the dimensions are the same
          m_matrix(:, :, nn, nkp) = mmn_tmp(:, :)
        end if
      end do
      deallocate (mmn_tmp, stat=ierr)
      if (ierr /= 0) call io_error('Error in deallocating mmn_tmp in overlap_read')
      close (mmn_in)
    endif

    if (disentanglement) then
!       call comms_bcast(m_matrix_orig(1,1,1,1),num_bands*num_bands*nntot*num_kpts)
      call comms_scatterv(m_matrix_orig_local, num_bands*num_bands*nntot*counts(my_node_id), &
                          m_matrix_orig, num_bands*num_bands*nntot*counts, num_bands*num_bands*nntot*displs)
    else
!       call comms_bcast(m_matrix(1,1,1,1),num_wann*num_wann*nntot*num_kpts)
      call comms_scatterv(m_matrix_local, num_wann*num_wann*nntot*counts(my_node_id), &
                          m_matrix, num_wann*num_wann*nntot*counts, num_wann*num_wann*nntot*displs)
    endif

    if (.not. use_bloch_phases) then
      if (on_root) then

        ! Read A_matrix from file wannier.amn
        amn_in = io_file_unit()
        open (unit=amn_in, file=trim(seedname)//'.amn', form='formatted', status='old', err=102)

        if (on_root) write (stdout, '(/a)', advance='no') ' Reading projections from '//trim(seedname)//'.amn : '

        ! Read the comment line
        read (amn_in, '(a)', err=104, end=104) dummy
        if (on_root) write (stdout, '(a)') trim(dummy)

        ! Read the number of bands, k-points and wannier functions
        read (amn_in, *, err=104, end=104) nb_tmp, nkp_tmp, np_tmp

        ! Checks
        if (nb_tmp .ne. num_bands) &
          call io_error(trim(seedname)//'.amn has not the right number of bands')
        if (nkp_tmp .ne. num_kpts) &
          call io_error(trim(seedname)//'.amn has not the right number of k-points')
        if (np_tmp .ne. num_proj) &
          call io_error(trim(seedname)//'.amn has not the right number of projections')

        if (num_proj > num_wann .and. .not. lselproj) &
          call io_error(trim(seedname)//'.amn has too many projections to be used without selecting a subset')

        ! Read the projections
        num_amn = num_bands*num_proj*num_kpts
        if (disentanglement) then
          do ncount = 1, num_amn
            read (amn_in, *, err=104, end=104) m, n, nkp, a_real, a_imag
            if (proj2wann_map(n) < 0) cycle
            a_matrix(m, proj2wann_map(n), nkp) = cmplx(a_real, a_imag, kind=dp)
          end do
        else
          do ncount = 1, num_amn
            read (amn_in, *, err=104, end=104) m, n, nkp, a_real, a_imag
            if (proj2wann_map(n) < 0) cycle
            u_matrix(m, proj2wann_map(n), nkp) = cmplx(a_real, a_imag, kind=dp)
          end do
        end if
        close (amn_in)
      endif

      if (disentanglement) then
        call comms_bcast(a_matrix(1, 1, 1), num_bands*num_wann*num_kpts)
      else
        call comms_bcast(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts)
      endif

    else

      do n = 1, num_kpts
        do m = 1, num_wann
          u_matrix(m, m, n) = cmplx_1
        end do
      end do

    end if

    ! If post-processing a Car-Parinello calculation (gamma only)
    ! then rotate M and A to the basis of Kohn-Sham eigenstates
    if (cp_pp) call overlap_rotate()

    ! Check Mmn(k,b) is symmetric in m and n for gamma_only case
!~      if (gamma_only) call overlap_check_m_symmetry()

    ! If we don't need to disentangle we can now convert from A to U
    ! And rotate M accordingly
![ysl-b]
!       if(.not.disentanglement .and. (.not.cp_pp) .and. (.not. use_bloch_phases )) &
!            call overlap_project
!~       if((.not.cp_pp) .and. (.not. use_bloch_phases )) then
!~         if (.not.disentanglement) then
!~            if ( .not. gamma_only ) then
!~               call overlap_project
!~            else
!~               call overlap_project_gamma()
!~            endif
!~         else
!~            if (gamma_only) call overlap_symmetrize()
!~         endif
!~       endif
!
!~[aam]
    if ((.not. disentanglement) .and. (.not. cp_pp) .and. (.not. use_bloch_phases)) then
      if (.not. gamma_only) then
        call overlap_project()
      else
        call overlap_project_gamma()
      endif
    endif
!~[aam]

    !~      if( gamma_only .and. use_bloch_phases ) then
    !~        write(stdout,'(1x,"+",76("-"),"+")')
    !~        write(stdout,'(3x,a)') 'WARNING: gamma_only and use_bloch_phases                 '
    !~        write(stdout,'(3x,a)') '         M must be calculated from *real* Bloch functions'
    !~        write(stdout,'(1x,"+",76("-"),"+")')
    !~      end if
![ysl-e]

    if (timing_level > 0) call io_stopwatch('overlap: read', 2)

    return
101 call io_error('Error: Problem opening input file '//trim(seedname)//'.mmn')
102 call io_error('Error: Problem opening input file '//trim(seedname)//'.amn')
103 call io_error('Error: Problem reading input file '//trim(seedname)//'.mmn')
104 call io_error('Error: Problem reading input file '//trim(seedname)//'.amn')

  end subroutine overlap_read

!~[aam]
!~  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!~ subroutine overlap_check_m_symmetry()
!~ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!~
!~   use w90_constants,  only : eps6
!~   use w90_parameters, only : num_bands,num_wann,a_matrix,m_matrix_orig,u_matrix,m_matrix, &
!~                              nntot,timing_level,disentanglement
!~   use w90_io,         only : io_error,io_stopwatch
!~
!~   implicit none
!~
!~   integer       :: nn,i,j,m,n, p(1), mdev, ndev, nndev,ierr, num_tmp
!~  real(kind=dp) :: dev,dev_tmp
!~
!~   if (timing_level>1) call io_stopwatch('overlap: check_m_sym',1)
!~
!~   if (disentanglement) then
!~      num_tmp=num_bands
!~   else
!~      num_tmp=num_wann
!~   endif
!~
!~   ! check whether M_mn is symmetric
!~   dev = 0.0_dp
!~   do nn=1,nntot
!~      do n=1,num_tmp
!~         do m=1,n
!~            if (disentanglement) then
!~               dev_tmp=abs(m_matrix_orig(m,n,nn,1)-m_matrix_orig(n,m,nn,1))
!~            else
!~               dev_tmp=abs(m_matrix(m,n,nn,1)-m_matrix(n,m,nn,1))
!~            endif
!~            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. eps6 ) then
!~      write(stdout,'(1x,"+",76("-"),"+")')
!~      write(stdout,'(3x,a)') 'WARNING: M is not strictly symmetric in overlap_check_m_symmetry'
!~      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
!~
!~   if (timing_level>1) call io_stopwatch('overlap: check_m_sym',2)
!~
!~   return
!~
!~ end subroutine overlap_check_m_symmetry
!~[aam]

!~![ysl-b]
!~  !%%%%%%%%%%%%%%%%%%%%%
!~  subroutine overlap_symmetrize
!~  !%%%%%%%%%%%%%%%%%%%%%
!~
!~    use w90_parameters, only : num_bands,num_wann,a_matrix,m_matrix_orig,u_matrix,m_matrix, &
!~                               nntot,timing_level,disentanglement
!~    use w90_io,         only : io_error,io_stopwatch
!~
!~    implicit none
!~
!~    integer       :: nn,i,j,m,n, p(1), mdev, ndev, nndev,ierr
!~    real(kind=dp) :: eps,dev,dev_tmp
!~    real(kind=dp),    allocatable :: a_cmp(:)
!~
!~
!~    if (timing_level>1) call io_stopwatch('overlap: symmetrize',1)
!~
!~    allocate(a_cmp(num_wann),stat=ierr)
!~    if (ierr/=0) call io_error('Error in allocating a_cmp in overlap_symmetrize')
!~    allocate(ph_g(num_bands),stat=ierr)
!~    if (ierr/=0) call io_error('Error in allocating ph_g in overlap_symmetrize')
!~
!~    eps = 1.0e-6_dp
!~    ph_g=cmplx_0
!~
!~    ! disentanglement - m_matrix_orig
!~    ! localzation only - m_matrix
!~
!~    ! If a wavefunction is real except for a phase factor e^(i*phi_m) = ph_g(m)
!~    ! A_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
!~    !
!~   do m=1,num_bands
!~       a_cmp(:)=abs(a_matrix(m,:,1))
!~       p=maxloc(a_cmp)
!~       ph_g(m)=conjg(a_matrix(m,p(1),1)/abs(a_matrix(m,p(1),1)))
!~   end do
!~
!~    ! M_mn (new) = ph_g(m) * M_mn (old) * conjg(ph_g(n))
!~    ! a_mn (new) = ph_g(m) * a_mn (old)
!~
!~    do nn=1,nntot
!~       do n=1,num_bands
!~          do m=1,num_bands
!~             m_matrix_orig(m,n,nn,1)=ph_g(m)*m_matrix_orig(m,n,nn,1)*conjg(ph_g(n))
!~          end do
!~       end do
!~    end do
!~    do n=1,num_wann
!~       do m=1,num_bands
!~          a_matrix(m,n,1)=ph_g(m)*a_matrix(m,n,1)
!~       end do
!~    end do
!~
!~    ! check whether M_mn is now symmetric
!~    dev = 0.0_dp
!~    do nn=1,nntot
!~       do n=1,num_bands
!~          do m=1,n
!~             dev_tmp=abs(m_matrix_orig(m,n,nn,1)-m_matrix_orig(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_symmetrize'
!~       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
!~
!~    deallocate(a_cmp,stat=ierr)
!~    if (ierr/=0) call io_error('Error in deallocating a_cmp in overlap_symmetrize')
!~
!~    if (timing_level>1) call io_stopwatch('overlap: symmetrize',2)
!~
!~    return
!~
!~  end subroutine overlap_symmetrize
!~![ysl-e]

  !%%%%%%%%%%%%%%%%%%%%%
  subroutine overlap_rotate
    !%%%%%%%%%%%%%%%%%%%%%
    !! Only used when interfaced to the CP code
    !! Not sure why this is done here and not in CP

    use w90_parameters, only: num_bands, a_matrix, m_matrix_orig, nntot, timing_level
    use w90_io, only: io_file_unit, io_error, io_stopwatch

    implicit none

    integer       :: lam_unit, info, inn, i, j
    real(kind=DP) :: lambda(num_bands, num_bands)
    real(kind=DP) :: AP(num_bands*(num_bands + 1)/2)
    real(kind=DP) :: eig(num_bands), work(3*num_bands)

    if (timing_level > 1) call io_stopwatch('overlap: rotate', 1)

    lam_unit = io_file_unit()
    open (unit=lam_unit, file='lambda.dat', &
          form='unformatted', status='old', action='read')
!~    write(stdout,*) ' Reading lambda.dat...'
    read (lam_unit) lambda
!~    write(stdout,*) ' done'
    close (lam_unit)

    do j = 1, num_bands
      do i = 1, j
        AP(i + (j - 1)*j/2) = 0.5_dp*(lambda(i, j) + lambda(j, i))
      end do
    end do

    CALL DSPEV('V', 'U', num_bands, AP, eig, lambda, num_bands, work, info)
    if (info .ne. 0) &
      call io_error('Diagonalization of lambda in overlap_rotate failed')

    ! For debugging
!~    write(stdout,*) 'EIGENVALUES - CHECK WITH CP OUTPUT'
!~    do i=1,num_bands
!~       write(stdout,*) 13.6058*eig(i)
!~    end do

    ! Rotate M_mn
    do inn = 1, nntot
      m_matrix_orig(:, :, inn, 1) = &
        matmul(transpose(lambda), matmul(m_matrix_orig(:, :, inn, 1), lambda))
    end do

    ! Rotate A_mn
    a_matrix(:, :, 1) = matmul(transpose(lambda), a_matrix(:, :, 1))

    ! For debugging
!~    ! Write rotated A and M
!~    do i=1,num_bands
!~       do j=1,num_wann
!~          write(12,'(2i5,a,2f18.12)') i,j,'   1',a_matrix(i,j,1)
!~       enddo
!~    enddo
!~    do inn=1,nntot
!~       do i=1,num_bands
!~          do j=1,num_wann
!~             write(11,'(2i5,2a)') i,j,'    1','    1'
!~             write(11,'(2f18.12)') m_matrix_orig(i,j,inn,1)
!~          enddo
!~       enddo
!~    enddo
!~    stop

    if (timing_level > 1) call io_stopwatch('overlap: rotate', 2)

    return

  end subroutine overlap_rotate

  !%%%%%%%%%%%%%%%%%%%%%
  subroutine overlap_dealloc()
    !%%%%%%%%%%%%%%%%%%%%%
    !! Dellocate memory

    use w90_parameters, only: u_matrix, m_matrix, m_matrix_orig, &
      a_matrix, u_matrix_opt, &
      m_matrix_local, m_matrix_orig_local
    use w90_io, only: io_error

    implicit none

    integer :: ierr

    if (allocated(u_matrix_opt)) then
      deallocate (u_matrix_opt, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating u_matrix_opt in overlap_dealloc')
    end if
    if (allocated(a_matrix)) then
      deallocate (a_matrix, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating a_matrix in overlap_dealloc')
    end if
    if (on_root) then
    if (allocated(m_matrix_orig)) then
      deallocate (m_matrix_orig, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating m_matrix_orig in overlap_dealloc')
    endif
    endif
    if (allocated(m_matrix_orig_local)) then
      deallocate (m_matrix_orig_local, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating m_matrix_orig_local in overlap_dealloc')
    endif
!~![ysl-b]
!~    if (allocated( ph_g)) then
!~       deallocate( ph_g, stat=ierr )
!~       if (ierr/=0) call io_error('Error deallocating ph_g in overlap_dealloc')
!~    endif
!~![ysl-e]

!    if (on_root) then
!    deallocate ( m_matrix, stat=ierr )
!    if (ierr/=0) call io_error('Error deallocating m_matrix in overlap_dealloc')
!    endif
!    deallocate ( m_matrix_local, stat=ierr )
!    if (ierr/=0) call io_error('Error deallocating m_matrix_local in overlap_dealloc')
!    deallocate ( u_matrix, stat=ierr )
!    if (ierr/=0) call io_error('Error deallocating u_matrix in overlap_dealloc')
    if (on_root) then
      if (allocated(m_matrix)) then
        deallocate (m_matrix, stat=ierr)
        if (ierr /= 0) call io_error('Error deallocating m_matrix in overlap_dealloc')
      endif
    endif
    if (allocated(m_matrix_local)) then
      deallocate (m_matrix_local, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating m_matrix_local in overlap_dealloc')
    endif
    if (allocated(u_matrix)) then
      deallocate (u_matrix, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating u_matrix in overlap_dealloc')
    endif

    return

  end subroutine overlap_dealloc

  !==================================================================!
  subroutine overlap_project()
    !==================================================================!
    !!  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_error, io_stopwatch
    use w90_parameters, only: num_bands, num_wann, num_kpts, timing_level, &
      u_matrix, m_matrix, nntot, nnlist, &
      m_matrix_local
    use w90_utility, only: utility_zgemm
    use w90_parameters, only: lsitesymmetry !RS:
    use w90_sitesym, only: sitesym_symmetrize_u_matrix !RS:
    use w90_comms, only: my_node_id, num_nodes, &
      comms_array_split, comms_scatterv, comms_gatherv

    implicit none

    ! internal variables
    integer :: i, j, m, nkp, 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(:, :)
    ! Needed to split an array on different nodes
    integer, dimension(0:num_nodes - 1) :: counts
    integer, dimension(0:num_nodes - 1) :: displs

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

    call comms_array_split(num_kpts, counts, displs)

    allocate (svals(num_bands), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating svals in overlap_project')
    allocate (cz(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating cz in overlap_project')
    allocate (cvdag(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating cvdag in overlap_project')
    allocate (cwork(4*num_bands), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating cwork in overlap_project')

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

    do nkp = 1, num_kpts
      !
      ! 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'
        endif
        call io_error('Error in ZGESVD in overlap_project')
      endif

!       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))
          enddo
          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 io_error('Error in unitarity of initial U in overlap_project')
          endif
          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 io_error('Error in unitarity of initial U in overlap_project')
          endif
        enddo
      enddo
    enddo
    ! NKP

    if (lsitesymmetry) call sitesym_symmetrize_u_matrix(num_wann, u_matrix) !RS: update U(Rk)

    ! 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 nkp = 1, counts(my_node_id)
      do nn = 1, nntot
        nkp2 = nnlist(nkp + displs(my_node_id), nn)
        ! cvdag = U^{dagger} . M   (use as workspace)
        call utility_zgemm(cvdag, u_matrix(:, :, nkp + displs(my_node_id)), 'C', m_matrix_local(:, :, nn, nkp), 'N', num_wann)
        ! cz = cvdag . U
        call utility_zgemm(cz, cvdag, 'N', u_matrix(:, :, nkp2), 'N', num_wann)
        m_matrix_local(:, :, nn, nkp) = cz(:, :)
      end do
    end do
    call comms_gatherv(m_matrix_local, num_wann*num_wann*nntot*counts(my_node_id), &
                       m_matrix, num_wann*num_wann*nntot*counts, num_wann*num_wann*nntot*displs)

    deallocate (cwork, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating cwork in overlap_project')
    deallocate (cvdag, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating cvdag in overlap_project')
    deallocate (cz, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating cz in overlap_project')
    deallocate (svals, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating svals in overlap_project')

    if (timing_level > 1) call io_stopwatch('overlap: project', 2)

    return

  end subroutine overlap_project

![ysl-b]
  !==================================================================!
  subroutine overlap_project_gamma()
    !==================================================================!
    !!  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_error, io_stopwatch
    use w90_parameters, only: num_wann, timing_level, &
      u_matrix, m_matrix, nntot!,num_kpts,nnlist
    use w90_utility, only: utility_zgemm

    implicit none

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

!~ real(kind=dp),    allocatable :: u_cmp(:)
!~ integer :: n,mdev, ndev, nndev,p(1)
!~ real(kind=dp)                 :: dev, dev_tmp

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

!~    allocate(ph_g(num_wann),stat=ierr)
!~    if (ierr/=0) call io_error('Error in allocating ph_g in overlap_project_gamma')
    ! internal variables
    allocate (u_matrix_r(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating u_matrix_r in overlap_project_gamma')
!~    allocate(u_cmp(num_wann),stat=ierr)
!~    if (ierr/=0) call io_error('Error in allocating u_cmp in overlap_project_gamma')
    allocate (svals(num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating svals in overlap_project_gamma')
    allocate (work(5*num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating work in overlap_project_gamma')
    allocate (rz(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating rz in overlap_project_gamma')
    allocate (rv(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating rv in overlap_project_gamma')
    allocate (cz(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating cz in overlap_project_gamma')
    allocate (cvdag(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating cvdag in overlap_project_gamma')

    !
!~    ! 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'
      endif
      call io_error('overlap_project_gamma: problem in DGESVD 1')
    endif

    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)
        enddo
        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 io_error('Error in unitarity of initial U in overlap_project_gamma')
        endif
        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 io_error('Error in unitarity of initial U in overlap_project_gamma')
        endif
      enddo
    enddo

    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) call io_error('Error in deallocating cvdag in overlap_project_gamma')
    deallocate (cz, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating cz in overlap_project_gamma')
    deallocate (rv, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating rv in overlap_project_gamma')
    deallocate (rz, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating rz in overlap_project_gamma')
    deallocate (work, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating work in overlap_project_gamma')
    deallocate (svals, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating svals in overlap_project_gamma')
!~    deallocate(u_cmp,stat=ierr)
!~    if (ierr/=0) call io_error('Error in deallocating u_cmp in overlap_project_gamma')
    deallocate (u_matrix_r, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating u_matrix_r in overlap_project_gamma')

    if (timing_level > 1) call io_stopwatch('overlap: project_gamma', 2)

    return

  end subroutine overlap_project_gamma
![ysl-e]

end module w90_overlap