Construct projections for the start of the disentanglement routine
Original notes from Nicola (refers only to the square case)
This subroutine calculates the transformation matrix CU = CS^(-1/2).CA, where CS = CA.CA^dagger. CS is diagonalized with a Schur factorization, to be on the safe side of numerical stability.
ZGEES computes for an N-by-N complex nonsymmetric matrix Y, the eigenvalues, the Schur form T, and, optionally, the matrix of Schur vectors. This gives the Schur factorization Y = ZT(Z**H).
Optionally, it also orders the eigenvalues on the diagonal of the Schur form so that selected eigenvalues are at the top left. The leading components of Z then form an orthonormal basis for the invariant subspace corresponding to the selected eigenvalues.
A complex matrix is in Schur form if it is upper triangular.
Notes from Ivo disentangling (e.g. non-square) projection (See Sec. III.D of SMV paper) Compute the ndimwin(k) x num_wann matrix cu that yields, from the ndimwin original Bloch states, the num_wann Bloch-like states with maximal projection onto the num_wann localised functions:
CU = CA.CS^{-1/2}, CS = transpose(CA).CA
Use the singular-calue decomposition of the matrix CA:
CA = CZ.CD.CVdagger (note: zgesvd spits out CVdagger)
which yields
CU = CZ.CD.CD^{-1}.CVdag
where CZ is ndimwin(NKP) x ndimwin(NKP) and unitary, CD is ndimwin(NKP) x num_wann and diagonal, CD^{-1} is num_wann x num_wann and diagonal, and CVdag is num_wann x num_wann and unitary.
subroutine dis_project()
!==================================================================!
! !
!! Construct projections for the start of the disentanglement routine
!!
!! Original notes from Nicola (refers only to the square case)
!!
!! This subroutine calculates the transformation matrix
!! CU = CS^(-1/2).CA, where CS = CA.CA^dagger.
!! CS is diagonalized with a Schur factorization, to be on the safe
!! side of numerical stability.
!!
!! ZGEES computes for an N-by-N complex nonsymmetric matrix Y, the
!! eigenvalues, the Schur form T, and, optionally, the matrix of
!! Schur vectors. This gives the Schur factorization Y = Z*T*(Z**H).
!!
!! Optionally, it also orders the eigenvalues on the diagonal of the
!! Schur form so that selected eigenvalues are at the top left.
!! The leading components of Z then form an orthonormal basis for
!! the invariant subspace corresponding to the selected eigenvalues.
!!
!! A complex matrix is in Schur form if it is upper triangular.
!!
!! Notes from Ivo disentangling (e.g. non-square) projection
!! (See Sec. III.D of SMV paper)
!! Compute the ndimwin(k) x num_wann matrix cu that yields,
!! from the ndimwin original Bloch states, the num_wann Bloch-like
!! states with maximal projection onto the num_wann localised
!! functions:
!!
!! CU = CA.CS^{-1/2}, CS = transpose(CA).CA
!!
!! Use the singular-calue decomposition of the matrix CA:
!!
!! CA = CZ.CD.CVdagger (note: zgesvd spits out CVdagger)
!!
!! which yields
!!
!! CU = CZ.CD.CD^{-1}.CVdag
!!
!! where CZ is ndimwin(NKP) x ndimwin(NKP) and unitary, CD is
!! ndimwin(NKP) x num_wann and diagonal, CD^{-1} is
!! num_wann x num_wann and diagonal, and CVdag is
!! num_wann x num_wann and unitary.
!!
!==================================================================!
use w90_constants, only: eps5
implicit none
! internal variables
integer :: i, j, l, m, nkp, info, ierr
real(kind=dp), allocatable :: svals(:)
real(kind=dp), allocatable :: rwork(:)
complex(kind=dp) :: ctmp2
complex(kind=dp), allocatable :: cwork(:)
complex(kind=dp), allocatable :: cz(:, :)
complex(kind=dp), allocatable :: cvdag(:, :)
! complex(kind=dp), allocatable :: catmpmat(:,:,:)
if (timing_level > 1) call io_stopwatch('dis: project', 1)
if (on_root) write (stdout, '(/1x,a)') &
' Unitarised projection of Wannier functions '
if (on_root) write (stdout, '(1x,a)') &
' ------------------------------------------ '
if (on_root) write (stdout, '(3x,a)') 'A_mn = <psi_m|g_n> --> S = A.A^+ --> U = S^-1/2.A'
if (on_root) write (stdout, '(3x,a)', advance='no') 'In dis_project...'
! allocate(catmpmat(num_bands,num_bands,num_kpts),stat=ierr)
! if (ierr/=0) call io_error('Error in allocating catmpmat in dis_project')
allocate (svals(num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating svals in dis_project')
allocate (rwork(5*num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rwork in dis_project')
allocate (cvdag(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cvdag in dis_project')
allocate (cz(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cz in dis_project')
allocate (cwork(4*num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cwork in dis_project')
! here we slim down the ca matrix
! up to here num_bands(=num_bands) X num_wann(=num_wann)
! do nkp = 1, num_kpts
! do j = 1, num_wann
! do i = 1, ndimwin(nkp)
! catmpmat(i,j,nkp) = a_matrix(nfirstwin(nkp)+i-1,j,nkp)
! enddo
! enddo
! do j = 1, num_wann
! a_matrix(1:ndimwin(nkp),j,nkp) = catmpmat(1:ndimwin(nkp),j,nkp)
! enddo
! do j = 1, num_wann
! a_matrix(ndimwin(nkp)+1:num_bands,j,nkp) = cmplx_0
! enddo
! enddo
! in order to reduce the memory usage, we don't use catmpmat.
do nkp = 1, num_kpts
if (ndimwin(nkp) .ne. num_bands) then
do j = 1, num_wann
do i = 1, ndimwin(nkp)
ctmp2 = a_matrix(nfirstwin(nkp) + i - 1, j, nkp)
a_matrix(i, j, nkp) = ctmp2
enddo
a_matrix(ndimwin(nkp) + 1:num_bands, j, nkp) = cmplx_0
enddo
endif
enddo
do nkp = 1, num_kpts
! SINGULAR VALUE DECOMPOSITION
call ZGESVD('A', 'A', ndimwin(nkp), num_wann, a_matrix(:, :, nkp), &
num_bands, svals, cz, num_bands, cvdag, num_bands, cwork, &
4*num_bands, rwork, info)
if (info .ne. 0) then
if (on_root) write (stdout, *) ' ERROR: IN ZGESVD IN dis_project'
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_project: problem in ZGESVD 1')
endif
! NOTE THAT - AT LEAST FOR LINUX MKL LAPACK - THE OUTPUT OF ZGESVD
! GIVES ALREADY Vdagger, SO A=Z.S.Vdagger IS ACTUALLY GIVEN BY cz.s.cvda
!
! also, cu is cz.cd.cd^-1.cvdag, and the asymmetric structure of cd.cd^-
! allows us to say cu=cz.cvdag, but where the sum on the inner index
! goes only from 1 to num_wann (i.e. DO l=1,num_wann). This is because
! cz.cd.cd^-1 is a matrix ndimwin(nkp) X num_wann, identical to the
! first num_wann columns of the ndimwin(nkp) X ndimwin(nkp) matrix cz
!
! same for ca: s is a ndimwin(nkp) X num_wann matrix that is
! zero everywhere but for its num_wann X num_wann top square part,
! that is diagonal. Multiplying cz by the s matrix is equivalent
! to moltiplying the first num_wann columns of cz, each by the correspondin
! diagonal element of s, that is s(L)
! I'm not sure why we reconstruct ca in what follows - in one explicit t
! [ aam: it is because a_matrix is overwritten by ZGESVD ]
! it seemed to be identical to the input ca (as it should be)
u_matrix_opt(:, :, nkp) = cmplx_0
a_matrix(:, :, nkp) = cmplx_0
do j = 1, num_wann
do i = 1, ndimwin(nkp)
do l = 1, num_wann
u_matrix_opt(i, j, nkp) = u_matrix_opt(i, j, nkp) + cz(i, l)*cvdag(l, j)
a_matrix(i, j, nkp) = a_matrix(i, j, nkp) + cz(i, l)*svals(l)*cvdag(l, j)
enddo
enddo
enddo
!
! CHECK UNITARITY
!
! note that cu.transpose(cu) is *NOT* an identity ndimwin(nkp) by ndimwi
! matrix, but transpose(cu).cu is a num_wann by num_wann identity matrix.
! I have once checked the former statement, now I will just leave here t
! for the latter (what this means is that the columns of cu are orthonor
! vectors).
do i = 1, num_wann
do j = 1, num_wann
ctmp2 = cmplx_0
do m = 1, ndimwin(nkp)
ctmp2 = ctmp2 + u_matrix_opt(m, j, nkp)*conjg(u_matrix_opt(m, i, nkp))
enddo
if ((i .eq. j) .and. (abs(ctmp2 - cmplx_1) .gt. eps5)) then
if (on_root) write (stdout, *) ' ERROR: unitarity of initial U'
if (on_root) write (stdout, '(1x,a,i2)') 'nkp= ', nkp
if (on_root) write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
if (on_root) write (stdout, '(1x,a,f12.6,1x,f12.6)') &
'[u_matrix_opt.transpose(u_matrix_opt)]_ij= ', &
real(ctmp2, dp), aimag(ctmp2)
call io_error('dis_project: Error in unitarity of initial U in dis_project')
endif
if ((i .ne. j) .and. (abs(ctmp2) .gt. eps5)) then
if (on_root) write (stdout, *) ' ERROR: unitarity of initial U'
if (on_root) write (stdout, '(1x,a,i2)') 'nkp= ', nkp
if (on_root) write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
if (on_root) write (stdout, '(1x,a,f12.6,1x,f12.6)') &
'[u_matrix_opt.transpose(u_matrix_opt)]_ij= ', &
real(ctmp2, dp), aimag(ctmp2)
call io_error('dis_project: Error in unitarity of initial U in dis_project')
endif
enddo
enddo
enddo
! NKP
deallocate (cwork, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cwork in dis_project')
deallocate (cz, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cz in dis_project')
deallocate (cvdag, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cvdag in dis_project')
deallocate (rwork, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating rwork in dis_project')
deallocate (svals, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating svals in dis_project')
! deallocate(catmpmat,stat=ierr)
! if (ierr/=0) call io_error('Error in deallocating catmpmat in dis_project')
if (on_root) write (stdout, '(a)') ' done'
if (timing_level > 1) call io_stopwatch('dis: project', 2)
return
end subroutine dis_project