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_read
  public :: overlap_dealloc
  public :: overlap_project
  public :: overlap_project_gamma  ![ysl]
!~  public :: overlap_check_m_symmetry

contains

    

  !%%%%%%%%%%%%%%%%%%%%%
  subroutine overlap_read( )
  !%%%%%%%%%%%%%%%%%%%%%
    !! Allocate and read the Mmn and Amn from files

    use w90_parameters, only : num_bands, num_wann, num_kpts, nntot, nncell, nnlist,&
                           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
    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, nw_tmp, ierr
    real(kind=dp) :: m_real, m_imag, a_real, a_imag
    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)

    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 (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, nw_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 (nw_tmp.ne.num_wann) &
               call io_error(trim(seedname)//'.amn has not the right number of Wannier functions')
          
          ! Read the projections
          num_amn = num_bands*num_wann*num_kpts
          if (disentanglement) then
             do ncount = 1, num_amn
                read(amn_in,*,err=104,end=104) m,n,nkp,a_real,a_imag
                a_matrix(m,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
                u_matrix(m,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