Main disentanglement routine
subroutine dis_main()
!==================================================================!
!! Main disentanglement routine
! !
! !
! !
!==================================================================!
use w90_io, only: io_file_unit
! internal variables
integer :: nkp, nkp2, nn, j, ierr, page_unit
integer :: nkp_global
complex(kind=dp), allocatable :: cwb(:, :), cww(:, :)
! 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('dis: main', 1)
call comms_array_split(num_kpts, counts, displs)
if (on_root) write (stdout, '(/1x,a)') &
'*------------------------------- DISENTANGLE --------------------------------*'
! Allocate arrays
allocate (eigval_opt(num_bands, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating eigval_opt in dis_main')
eigval_opt = eigval
! Set up energy windows
call dis_windows()
! Construct the unitarized projection
call dis_project()
! If there is an inner window, need to modify projection procedure
! (Sec. III.G SMV)
if (linner) then
if (lsitesymmetry) call io_error('in symmetry-adapted mode, frozen window not implemented yet') !YN: RS:
if (on_root) write (stdout, '(3x,a)') 'Using an inner window (linner = T)'
call dis_proj_froz()
else
if (on_root) write (stdout, '(3x,a)') 'No inner window (linner = F)'
endif
! Debug
call internal_check_orthonorm()
! Slim down the original Mmn(k,b)
call internal_slim_m()
lwindow = .false.
do nkp = 1, num_kpts
do j = nfirstwin(nkp), nfirstwin(nkp) + ndimwin(nkp) - 1
lwindow(j, nkp) = .true.
end do
end do
if (lsitesymmetry) call sitesym_slim_d_matrix_band(lwindow) !RS: calculate initial U_{opt}(Rk) from U_{opt}(k)
if (lsitesymmetry) call sitesym_symmetrize_u_matrix(num_bands, u_matrix_opt, lwindow) !RS:
! Extract the optimally-connected num_wann-dimensional subspaces
![ysl-b]
if (.not. gamma_only) then
call dis_extract()
else
call dis_extract_gamma()
end if
![ysl-e]
! Allocate workspace
allocate (cwb(num_wann, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cwb in dis_main')
allocate (cww(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cww in dis_main')
! Find the num_wann x num_wann overlap matrices between
! the basis states of the optimal subspaces
do nkp = 1, counts(my_node_id)
nkp_global = nkp + displs(my_node_id)
do nn = 1, nntot
nkp2 = nnlist(nkp_global, nn)
call zgemm('C', 'N', num_wann, ndimwin(nkp2), ndimwin(nkp_global), cmplx_1, &
u_matrix_opt(:, :, nkp_global), num_bands, m_matrix_orig_local(:, :, nn, nkp), num_bands, &
cmplx_0, cwb, num_wann)
call zgemm('N', 'N', num_wann, num_wann, ndimwin(nkp2), cmplx_1, &
cwb, num_wann, u_matrix_opt(:, :, nkp2), num_bands, &
cmplx_0, cww, num_wann)
m_matrix_orig_local(1:num_wann, 1:num_wann, nn, nkp) = cww(:, :)
enddo
enddo
! Find the initial u_matrix
if (lsitesymmetry) call sitesym_replace_d_matrix_band() !RS: replace d_matrix_band here
![ysl-b]
if (.not. gamma_only) then
call internal_find_u()
else
call internal_find_u_gamma()
end if
![ysl-e]
if (optimisation <= 0) then
page_unit = io_file_unit()
open (unit=page_unit, form='unformatted', status='scratch')
! Update the m_matrix accordingly
do nkp = 1, counts(my_node_id)
nkp_global = nkp + displs(my_node_id)
do nn = 1, nntot
nkp2 = nnlist(nkp_global, nn)
call zgemm('C', 'N', num_wann, num_wann, num_wann, cmplx_1, &
u_matrix(:, :, nkp_global), num_wann, m_matrix_orig_local(:, :, nn, nkp), num_bands, &
cmplx_0, cwb, num_wann)
call zgemm('N', 'N', num_wann, num_wann, num_wann, cmplx_1, &
cwb, num_wann, u_matrix(:, :, nkp2), num_wann, &
cmplx_0, cww, num_wann)
write (page_unit) cww(:, :)
enddo
enddo
rewind (page_unit)
deallocate (m_matrix_orig_local, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating m_matrix_orig_local in dis_main')
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 dis_main')
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 dis_main')
do nkp = 1, counts(my_node_id)
do nn = 1, nntot
read (page_unit) m_matrix_local(:, :, nn, nkp)
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)
close (page_unit)
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 dis_main')
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 dis_main')
! Update the m_matrix accordingly
do nkp = 1, counts(my_node_id)
nkp_global = nkp + displs(my_node_id)
do nn = 1, nntot
nkp2 = nnlist(nkp_global, nn)
call zgemm('C', 'N', num_wann, num_wann, num_wann, cmplx_1, &
u_matrix(:, :, nkp_global), num_wann, m_matrix_orig_local(:, :, nn, nkp), num_bands, &
cmplx_0, cwb, num_wann)
call zgemm('N', 'N', num_wann, num_wann, num_wann, cmplx_1, &
cwb, num_wann, u_matrix(:, :, nkp2), num_wann, &
cmplx_0, cww, num_wann)
m_matrix_local(:, :, nn, nkp) = cww(:, :)
enddo
enddo
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 (m_matrix_orig_local, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating m_matrix_orig_local in dis_main')
endif
deallocate (a_matrix, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating a_matrix in dis_main')
! Deallocate workspace
deallocate (cww, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cww in dis_main')
deallocate (cwb, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cwb in dis_main')
!zero the unused elements of u_matrix_opt (just in case...)
do nkp = 1, num_kpts
do j = 1, num_wann
if (ndimwin(nkp) < num_bands) &
u_matrix_opt(ndimwin(nkp) + 1:, j, nkp) = cmplx_0
enddo
enddo
!~![ysl-b]
!~! Apply phase factor ph_g if gamma_only
!~ if (.not. gamma_only) then
!~ do nkp = 1, num_kpts
!~ do j = 1, num_wann
!~ u_matrix_opt(1:ndimwin(nkp),j,nkp) = u_matrix_opt(1:ndimwin(nkp),j,nkp)
!~ enddo
!~ enddo
!~ else
!~ do nkp = 1, num_kpts
!~ do j = 1, ndimwin(nkp)
!~ u_matrix_opt(j,1:num_wann,nkp) = conjg(ph_g(j))*u_matrix_opt(j,1:num_wann,nkp)
!~ enddo
!~ enddo
!~ endif
!~![ysl-e]
! Deallocate module arrays
call internal_dealloc()
if (timing_level > 0 .and. on_root) call io_stopwatch('dis: main', 2)
return
contains
!================================================================!
subroutine internal_check_orthonorm()
!================================================================!
! !
!! This subroutine checks that the states in the columns of the
!! final matrix U_opt are orthonormal at every k-point, i.e.,
!! that the matrix is unitary in the sense that
!! conjg(U_opt).U_opt = 1 (but not U_opt.conjg(U_opt) = 1).
!!
!! In particular, this checks whether the projected gaussians
!! are indeed orthogonal to the frozen states, at those k-points
!! where both are present in the trial subspace.
! !
!================================================================!
use w90_constants, only: eps8
implicit none
integer :: nkp, l, m, j
complex(kind=dp) :: ctmp
if (timing_level > 1) call io_stopwatch('dis: main: check_orthonorm', 1)
do nkp = 1, num_kpts
do l = 1, num_wann
do m = 1, l
ctmp = cmplx_0
do j = 1, ndimwin(nkp)
ctmp = ctmp + conjg(u_matrix_opt(j, m, nkp))*u_matrix_opt(j, l, nkp)
enddo
if (l .eq. m) then
if (abs(ctmp - cmplx_1) .gt. eps8) then
if (on_root) write (stdout, '(3i6,2f16.12)') nkp, l, m, ctmp
if (on_root) write (stdout, '(1x,a)') 'The trial orbitals for disentanglement are not orthonormal'
! write(stdout,'(1x,a)') 'Try re-running the calculation with the input keyword'
! write(stdout,'(1x,a)') ' devel_flag=orth-fix'
! write(stdout,'(1x,a)') 'Please report the sucess or failure of this to the Wannier90 developers'
call io_error('Error in dis_main: orthonormal error 1')
endif
else
if (abs(ctmp) .gt. eps8) then
if (on_root) write (stdout, '(3i6,2f16.12)') nkp, l, m, ctmp
if (on_root) write (stdout, '(1x,a)') 'The trial orbitals for disentanglement are not orthonormal'
! write(stdout,'(1x,a)') 'Try re-running the calculation with the input keyword'
! write(stdout,'(1x,a)') ' devel_flag=orth-fix'
! write(stdout,'(1x,a)') 'Please report the sucess or failure of this to the Wannier90 developers'
call io_error('Error in dis_main: orthonormal error 2')
endif
endif
enddo
enddo
enddo
if (timing_level > 1 .and. on_root) call io_stopwatch('dis: main: check_orthonorm', 2)
return
end subroutine internal_check_orthonorm
!================================================================!
subroutine internal_slim_m()
!================================================================!
! !
!! This subroutine slims down the original Mmn(k,b), removing
!! rows and columns corresponding to u_nks that fall outside
!! the outer energy window.
! !
!================================================================!
implicit none
integer :: nkp, nkp2, nn, i, j, m, n, ierr
integer :: nkp_global
complex(kind=dp), allocatable :: cmtmp(:, :)
! 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 .and. on_root) call io_stopwatch('dis: main: slim_m', 1)
call comms_array_split(num_kpts, counts, displs)
allocate (cmtmp(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cmtmp in dis_main')
do nkp = 1, counts(my_node_id)
nkp_global = nkp + displs(my_node_id)
do nn = 1, nntot
nkp2 = nnlist(nkp_global, nn)
do j = 1, ndimwin(nkp2)
n = nfirstwin(nkp2) + j - 1
do i = 1, ndimwin(nkp_global)
m = nfirstwin(nkp_global) + i - 1
cmtmp(i, j) = m_matrix_orig_local(m, n, nn, nkp)
enddo
enddo
m_matrix_orig_local(:, :, nn, nkp) = cmplx_0
do j = 1, ndimwin(nkp2)
do i = 1, ndimwin(nkp_global)
m_matrix_orig_local(i, j, nn, nkp) = cmtmp(i, j)
enddo
enddo
enddo
enddo
deallocate (cmtmp, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cmtmp in dis_main')
if (timing_level > 1 .and. on_root) call io_stopwatch('dis: main: slim_m', 2)
return
end subroutine internal_slim_m
!================================================================!
subroutine internal_find_u()
!================================================================!
! !
!! This subroutine finds the initial guess for the square unitary
!! rotation matrix u_matrix. The method is similar to Sec. III.D
!! of SMV, but with square instead of rectangular matrices:
!!
!! First find caa, the square overlap matrix <psitilde_nk|g_m>,
!! where psitilde is an eigenstate of the optimal subspace.
!!
!! Note that, contrary to what is implied in Sec. III.E of SMV,
!! this does *not* need to be computed by brute: instead we take
!! advantage of the previous computation of overlaps with the
!! same projections that are used to initiate the minimization of
!! Omega.
!!
!! Note: |psi> U_opt = |psitilde> and obviously
!! <psitilde| = (U_opt)^dagger <psi|
! !
!================================================================!
use w90_sitesym, only: ir2ik, ik2ir !YN: RS:
implicit none
integer :: nkp, info, ierr
complex(kind=dp), allocatable :: caa(:, :, :)
! For ZGESVD
real(kind=dp), allocatable :: svals(:)
real(kind=dp), allocatable :: rwork(:)
complex(kind=dp), allocatable :: cv(:, :)
complex(kind=dp), allocatable :: cz(:, :)
complex(kind=dp), allocatable :: cwork(:)
if (timing_level > 1 .and. on_root) call io_stopwatch('dis: main: find_u', 1)
! Currently, this part is not parallelized; thus, we perform the task only on root and then broadcast the result.
if (on_root) then
! Allocate arrays needed for ZGESVD
allocate (svals(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating svals in dis_main')
allocate (rwork(5*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rwork in dis_main')
allocate (cv(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cv in dis_main')
allocate (cz(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cz in dis_main')
allocate (cwork(4*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cwork in dis_main')
allocate (caa(num_wann, num_wann, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating caa in dis_main')
do nkp = 1, num_kpts
if (lsitesymmetry) then !YN: RS:
if (ir2ik(ik2ir(nkp)) .ne. nkp) cycle !YN: RS:
endif !YN: RS:
call zgemm('C', 'N', num_wann, num_wann, ndimwin(nkp), cmplx_1, &
u_matrix_opt(:, :, nkp), num_bands, a_matrix(:, :, nkp), num_bands, &
cmplx_0, caa(:, :, nkp), num_wann)
! Singular-value decomposition
call ZGESVD('A', 'A', num_wann, num_wann, caa(:, :, nkp), num_wann, &
svals, cz, num_wann, cv, num_wann, cwork, 4*num_wann, rwork, info)
if (info .ne. 0) then
if (on_root) write (stdout, *) ' ERROR: IN ZGESVD IN dis_main'
if (on_root) write (stdout, *) 'K-POINT NKP=', nkp, ' INFO=', info
if (info .lt. 0) then
if (on_root) write (stdout, *) 'THE ', -info, '-TH ARGUMENT HAD ILLEGAL VALUE'
endif
call io_error('dis_main: problem in ZGESVD 1')
endif
! u_matrix is the initial guess for the unitary rotation of the
! basis states given by the subroutine extract
call zgemm('N', 'N', num_wann, num_wann, num_wann, cmplx_1, &
cz, num_wann, cv, num_wann, cmplx_0, u_matrix(:, :, nkp), num_wann)
enddo
endif
call comms_bcast(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts)
! if (lsitesymmetry) call sitesym_symmetrize_u_matrix(num_wann,u_matrix) !RS:
if (on_root) then
! Deallocate arrays for ZGESVD
deallocate (caa, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating caa in dis_main')
deallocate (cwork, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cwork in dis_main')
deallocate (cz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cz in dis_main')
deallocate (cv, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cv in dis_main')
deallocate (rwork, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating rwork in dis_main')
deallocate (svals, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating svals in dis_main')
endif
if (lsitesymmetry) call sitesym_symmetrize_u_matrix(num_wann, u_matrix) !RS:
if (timing_level > 1) call io_stopwatch('dis: main: find_u', 2)
return
end subroutine internal_find_u
![ysl-b]
!================================================================!
subroutine internal_find_u_gamma()
!================================================================!
! !
!! Make initial u_matrix real
!! Must be the case when gamma_only = .true.
! !
!================================================================!
implicit none
integer :: info, ierr
real(kind=dp), allocatable :: u_opt_r(:, :)
real(kind=dp), allocatable :: a_matrix_r(:, :)
real(kind=dp), allocatable :: raa(:, :)
! For DGESVD
real(kind=dp), allocatable :: svals(:)
real(kind=dp), allocatable :: work(:)
real(kind=dp), allocatable :: rv(:, :)
real(kind=dp), allocatable :: rz(:, :)
if (timing_level > 1) call io_stopwatch('dis: main: find_u_gamma', 1)
! Allocate arrays needed for getting a_matrix_r
allocate (u_opt_r(ndimwin(1), num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating u_opt_r in dis_main')
allocate (a_matrix_r(ndimwin(1), num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating a_matrix_r in dis_main')
! Allocate arrays needed for DGESVD
allocate (svals(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating svals in dis_main')
allocate (work(5*num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rwork in dis_main')
allocate (rv(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cv in dis_main')
allocate (rz(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cz in dis_main')
allocate (raa(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating raa in dis_main')
u_opt_r(:, :) = real(u_matrix_opt(1:ndimwin(1), 1:num_wann, 1), dp)
a_matrix_r(:, :) = real(a_matrix(1:ndimwin(1), 1:num_wann, 1), kind=dp)
call dgemm('T', 'N', num_wann, num_wann, ndimwin(1), 1.0_dp, &
u_opt_r, ndimwin(1), a_matrix_r, ndimwin(1), &
0.0_dp, raa, num_wann)
! Singular-value decomposition
call DGESVD('A', 'A', num_wann, num_wann, raa, 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 dis_main'
write (stdout, *) 'K-POINT = Gamma', ' INFO=', info
if (info .lt. 0) then
write (stdout, *) 'THE ', -info, '-TH ARGUMENT HAD ILLEGAL VALUE'
endif
call io_error('dis_main: problem in DGESVD 1')
endif
! u_matrix is the initial guess for the unitary rotation of the
! basis states given by the subroutine extract
call dgemm('N', 'N', num_wann, num_wann, num_wann, 1.0_dp, &
rz, num_wann, rv, num_wann, 0.0_dp, raa, num_wann)
u_matrix(:, :, 1) = cmplx(raa(:, :), 0.0_dp, dp)
! Deallocate arrays for DGESVD
deallocate (raa, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating raa in dis_main')
deallocate (rz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating rz in dis_main')
deallocate (rv, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating rv in dis_main')
deallocate (work, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating work in dis_main')
deallocate (svals, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating svals in dis_main')
! Deallocate arrays for a_matrix_r
deallocate (a_matrix_r, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating a_matrix_r in dis_main')
deallocate (u_opt_r, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating u_opt_r in dis_main')
if (timing_level > 1) call io_stopwatch('dis: main: find_u_gamma', 2)
return
end subroutine internal_find_u_gamma
![ysl-e]
!==================================!
subroutine internal_dealloc()
!==================================!
!! Deallocate module data
! !
!==================================!
implicit none
integer :: ierr
! Module arrays allocated in dis_windows
deallocate (lfrozen, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating lfrozen in dis_main')
deallocate (indxnfroz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating indxnfroz in dis_main')
deallocate (indxfroz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating indxfroz in dis_main')
deallocate (ndimfroz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating ndimfroz in dis_main')
deallocate (nfirstwin, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating nfirstwin in dis_main')
! Module arrays allocated in dis_main
deallocate (eigval_opt, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating eigval_opt in dis_main')
return
end subroutine internal_dealloc
end subroutine dis_main