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
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