subroutine wann_svd_omega_i()
!========================================!
use w90_constants, only: dp, cmplx_0
use w90_io, only: io_stopwatch, io_error, stdout
use w90_parameters, only: num_wann, num_kpts, nntot, wb, &
m_matrix, lenconfac, length_unit, &
timing_level
implicit none
complex(kind=dp), allocatable :: cv1(:, :), cv2(:, :)
complex(kind=dp), allocatable :: cw1(:), cw2(:)
complex(kind=dp), allocatable :: cpad1(:)
real(kind=dp), allocatable :: singvd(:)
integer :: ierr, info
integer :: nkp, nn, nb, na, ind
real(kind=dp) :: omt1, omt2, omt3
if (timing_level > 1 .and. on_root) call io_stopwatch('wann: svd_omega_i', 1)
allocate (cw1(10*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cw1 in wann_svd_omega_i')
allocate (cw2(10*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cw2 in wann_svd_omega_i')
allocate (cv1(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cv1 in wann_svd_omega_i')
allocate (cv2(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cv2 in wann_svd_omega_i')
allocate (singvd(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating singvd in wann_svd_omega_i')
allocate (cpad1(num_wann*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cpad1 in wann_svd_omega_i')
cw1 = cmplx_0; cw2 = cmplx_0; cv1 = cmplx_0; cv2 = cmplx_0; cpad1 = cmplx_0
singvd = 0.0_dp
! singular value decomposition
omt1 = 0.0_dp; omt2 = 0.0_dp; omt3 = 0.0_dp
do nkp = 1, num_kpts
do nn = 1, nntot
ind = 1
do nb = 1, num_wann
do na = 1, num_wann
cpad1(ind) = m_matrix(na, nb, nn, nkp)
ind = ind + 1
enddo
enddo
call zgesvd('A', 'A', num_wann, num_wann, cpad1, num_wann, singvd, cv1, &
num_wann, cv2, num_wann, cw1, 10*num_wann, cw2, info)
if (info .ne. 0) then
call io_error('ERROR: Singular value decomp. zgesvd failed')
endif
do nb = 1, num_wann
omt1 = omt1 + wb(nn)*(1.0_dp - singvd(nb)**2)
omt2 = omt2 - wb(nn)*(2.0_dp*log(singvd(nb)))
omt3 = omt3 + wb(nn)*(acos(singvd(nb))**2)
enddo
enddo
enddo
omt1 = omt1/real(num_kpts, dp)
omt2 = omt2/real(num_kpts, dp)
omt3 = omt3/real(num_kpts, dp)
if (on_root) then
write (stdout, *) ' '
write (stdout, '(2x,a,f15.9,1x,a)') 'Omega Invariant: 1-s^2 = ', &
omt1*lenconfac**2, '('//trim(length_unit)//'^2)'
write (stdout, '(2x,a,f15.9,1x,a)') ' -2log s = ', &
omt2*lenconfac**2, '('//trim(length_unit)//'^2)'
write (stdout, '(2x,a,f15.9,1x,a)') ' acos^2 = ', &
omt3*lenconfac**2, '('//trim(length_unit)//'^2)'
endif
deallocate (cpad1, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cpad1 in wann_svd_omega_i')
deallocate (singvd, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating singvd in wann_svd_omega_i')
deallocate (cv2, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cv2 in wann_svd_omega_i')
deallocate (cv1, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cv1 in wann_svd_omega_i')
deallocate (cw2, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cw2 in wann_svd_omega_i')
deallocate (cw1, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cw1 in wann_svd_omega_i')
if (timing_level > 1 .and. on_root) call io_stopwatch('wann: svd_omega_i', 2)
return
end subroutine wann_svd_omega_i