subroutine plot_svd_omega_i(num_wann, num_kpts, kmesh_info, m_matrix, print_output, timer, &
dist_k, error, comm, stdout)
!================================================!
use w90_comms, only: w90_comm_type, mpirank, comms_allreduce
use w90_constants, only: dp, cmplx_0
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_types, only: kmesh_info_type, print_output_type, timer_list_type
implicit none
! arguments
complex(kind=dp), intent(in) :: m_matrix(:, :, :, :)
integer, intent(in) :: dist_k(:)
integer, intent(in) :: num_wann, num_kpts
integer, intent(in) :: stdout
type(kmesh_info_type), intent(in) :: kmesh_info
type(print_output_type), intent(in) :: print_output
type(timer_list_type), intent(inout) :: timer
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
! local variables
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, nkrank
integer :: nkp, nn, nb, na, ind
real(kind=dp) :: omt1, omt2, omt3
nkrank = count(dist_k == mpirank(comm)) ! number k this rank, for dimensioning
if (print_output%timing_level > 1 .and. print_output%iprint > 0) then
call io_stopwatch_start('wann: svd_omega_i', timer)
end if
allocate (cw1(10*num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cw1 in plot_svd_omega_i', comm)
return
end if
allocate (cw2(10*num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cw2 in plot_svd_omega_i', comm)
return
end if
allocate (cv1(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cv1 in plot_svd_omega_i', comm)
return
end if
allocate (cv2(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cv2 in plot_svd_omega_i', comm)
return
end if
allocate (singvd(num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating singvd in plot_svd_omega_i', comm)
return
end if
allocate (cpad1(num_wann*num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cpad1 in plot_svd_omega_i', comm)
return
end if
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, nkrank
do nn = 1, kmesh_info%nntot
ind = 1
do nb = 1, num_wann
do na = 1, num_wann
cpad1(ind) = m_matrix(na, nb, nn, nkp)
ind = ind + 1
end do
end do
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 set_error_fatal(error, 'ERROR: Singular value decomp. zgesvd failed', comm)
return
end if
do nb = 1, num_wann
omt1 = omt1 + kmesh_info%wb(nn)*(1.0_dp - singvd(nb)**2)
omt2 = omt2 - kmesh_info%wb(nn)*(2.0_dp*log(singvd(nb)))
omt3 = omt3 + kmesh_info%wb(nn)*(acos(singvd(nb))**2)
end do
end do
end do
call comms_allreduce(omt1, 1, 'SUM', error, comm)
call comms_allreduce(omt2, 1, 'SUM', error, comm)
call comms_allreduce(omt3, 1, 'SUM', error, comm)
omt1 = omt1/real(num_kpts, dp)
omt2 = omt2/real(num_kpts, dp)
omt3 = omt3/real(num_kpts, dp)
if (print_output%iprint > 0) then
write (stdout, *) ' '
write (stdout, '(2x,a,f15.9,1x,a)') 'Omega Invariant: 1-s^2 = ', &
omt1*print_output%lenconfac**2, '('//trim(print_output%length_unit)//'^2)'
write (stdout, '(2x,a,f15.9,1x,a)') ' -2log s = ', &
omt2*print_output%lenconfac**2, '('//trim(print_output%length_unit)//'^2)'
write (stdout, '(2x,a,f15.9,1x,a)') ' acos^2 = ', &
omt3*print_output%lenconfac**2, '('//trim(print_output%length_unit)//'^2)'
end if
deallocate (cpad1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cpad1 in plot_svd_omega_i', comm)
return
end if
deallocate (singvd, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating singvd in plot_svd_omega_i', comm)
return
end if
deallocate (cv2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cv2 in plot_svd_omega_i', comm)
return
end if
deallocate (cv1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cv1 in plot_svd_omega_i', comm)
return
end if
deallocate (cw2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cw2 in plot_svd_omega_i', comm)
return
end if
deallocate (cw1, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cw1 in plot_svd_omega_i', comm)
return
end if
if (print_output%timing_level > 1 .and. print_output%iprint > 0) then
call io_stopwatch_stop('wann: svd_omega_i', timer)
end if
end subroutine plot_svd_omega_i
