COMPUTES THE LEADING EIGENVECTORS OF Q_froz . P_s . Q_froz, WHERE P_s PROJECTOR OPERATOR ONTO THE SUBSPACE S OF THE PROJECTED GAUSSIANS, P_f THE PROJECTOR ONTO THE FROZEN STATES, AND Q_froz = 1 - P_froz, ALL EXP IN THE BASIS OF THE BLOCH EIGENSTATES INSIDE THE OUTER ENERGY WINDOW (See Eq. (27) in Sec. III.G of SMV)
subroutine dis_proj_froz()
!==================================================================!
! !
!! COMPUTES THE LEADING EIGENVECTORS OF Q_froz . P_s . Q_froz,
!! WHERE P_s PROJECTOR OPERATOR ONTO THE SUBSPACE S OF THE PROJECTED
!! GAUSSIANS, P_f THE PROJECTOR ONTO THE FROZEN STATES, AND
!! Q_froz = 1 - P_froz, ALL EXP IN THE BASIS OF THE BLOCH
!! EIGENSTATES INSIDE THE OUTER ENERGY WINDOW
!! (See Eq. (27) in Sec. III.G of SMV)
! !
!==================================================================!
use w90_constants, only: eps8
implicit none
! INPUT: num_wann,ndimwin,ndimfroz,indxfroz,lfrozen
! MODIFIED: u_matrix_opt (At input it contains the gaussians projected onto
! the window states in the routine project.f. At output
! the entries with the second index from 1 to ndimfroz(nkp)
! contain the frozen (inner window) states, while those
! from ndimfroz(nkp)+1 to num_wann have been replaced by
! the new trial states outside the inner window.)
! *********************************************************
! VARIABLES USED BY LAPACK'S ZHPEVX DIAGONALIZATION ROUTINE
! *********************************************************
integer, allocatable :: iwork(:)
integer, allocatable :: ifail(:)
real(kind=dp), allocatable :: w(:)
real(kind=dp), allocatable :: rwork(:)
complex(kind=dp), allocatable :: cap(:)
complex(kind=dp), allocatable :: cwork(:)
complex(kind=dp), allocatable :: cz(:, :)
! *********
! INTERNAL:
! *********
!
! CP_S(M,N) PROJECTION OPERATOR ONTO THE SUBSPACE OF THE PROJEC
! GAUSSIANS, EXPRESSED IN THE BASIS OF THE ORIGINAL BL
! EIGENSTATES INSIDE THE OUTER WINDOW (FOR THE PRESENT
! K-POINT)
! CQ_FROZ(M,N) PROJECTION OPERATOR ONTO THE SUBSPACE OF THE STATES
! THE SPACE OF FROZEN STATES (BUT INSIDE THE OUTER WIN
! EXPRESSED IN THE BASIS OF THE ORIGINAL BLOCH EIGENST
! INSIDE THE OUTER WINDOW (FOR THE PRESENT K-POINT)
! CPQ(M,N) THE MATRIX cp_s . cq_froz FOR THE PRESENT K-POINT
! CQPQ(M,N) THE MATRIX cq_froz . cp_s . cq_froz FOR THE PRESENT
!
integer :: goods, il, iu, nkp, l, j, n, m, info, ierr
integer :: counter, loop_f, loop_v, vmap(num_bands)
integer :: nzero
logical :: take
character(len=4) :: rep
complex(kind=dp) :: ctmp
complex(kind=dp), allocatable :: cp_s(:, :)
complex(kind=dp), allocatable :: cq_froz(:, :)
complex(kind=dp), allocatable :: cpq(:, :)
complex(kind=dp), allocatable :: cqpq(:, :)
if (timing_level > 1) call io_stopwatch('dis: proj_froz', 1)
if (on_root) write (stdout, '(3x,a)', advance='no') 'In dis_proj_froz...'
allocate (iwork(5*num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating iwork in dis_proj_froz')
allocate (ifail(num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating ifail in dis_proj_froz')
allocate (w(num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating w in dis_proj_froz')
allocate (rwork(7*num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating rwork in dis_proj_froz')
allocate (cap((num_bands*(num_bands + 1))/2), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cap in dis_proj_froz')
allocate (cwork(2*num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cwork in dis_proj_froz')
allocate (cz(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cz in dis_proj_froz')
allocate (cp_s(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cp_s in dis_proj_froz')
allocate (cq_froz(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cq_froz in dis_proj_froz')
allocate (cpq(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cpq in dis_proj_froz')
allocate (cqpq(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error allocating cqpq in dis_proj_froz')
do nkp = 1, num_kpts
! aam: this should be done at the end, otherwise important
! projection info is lost
!~ ! Put the frozen states in the lowest columns of u_matrix_opt
!~ if (ndimfroz(nkp).gt.0) then
!~ do l = 1, ndimfroz(nkp)
!~ u_matrix_opt(:,l,nkp)=cmplx_0
!~ u_matrix_opt(indxfroz(l,nkp),l,nkp) = cmplx_1
!~ enddo
!~ endif
! If there are non-frozen states, compute the num_wann-ndimfroz(nkp) leadin
! eigenvectors of cqpq
if (num_wann .gt. ndimfroz(nkp)) then
cq_froz = cmplx_0
cp_s = cmplx_0
do n = 1, ndimwin(nkp)
do m = 1, ndimwin(nkp)
do l = 1, num_wann
cp_s(m, n) = cp_s(m, n) + u_matrix_opt(m, l, nkp)*conjg(u_matrix_opt(n, l, nkp))
enddo
enddo
if (.not. lfrozen(n, nkp)) cq_froz(n, n) = cmplx_1
enddo
cpq = cmplx_0
do n = 1, ndimwin(nkp)
do m = 1, ndimwin(nkp)
do l = 1, ndimwin(nkp)
cpq(m, n) = cpq(m, n) + cp_s(m, l)*cq_froz(l, n)
enddo
enddo
enddo
cqpq = cmplx_0
do n = 1, ndimwin(nkp)
do m = 1, ndimwin(nkp)
do l = 1, ndimwin(nkp)
cqpq(m, n) = cqpq(m, n) + cq_froz(m, l)*cpq(l, n)
enddo
enddo
enddo
! DEBUG
! check hermiticity of cqpq
do n = 1, ndimwin(nkp)
do m = 1, n
if (abs(cqpq(m, n) - conjg(cqpq(n, m))) .gt. eps8) then
if (on_root) write (stdout, *) ' matrix CQPQ is not hermitian'
if (on_root) write (stdout, *) ' k-point ', nkp
call io_error('dis_proj_froz: error')
endif
enddo
enddo
! ENDDEBUG
cap = cmplx_0
do n = 1, ndimwin(nkp)
do m = 1, n
cap(m + (n - 1)*n/2) = cqpq(m, n)
enddo
enddo
il = ndimwin(nkp) - (num_wann - ndimfroz(nkp)) + 1
iu = ndimwin(nkp)
call ZHPEVX('V', 'A', 'U', ndimwin(nkp), cap, 0.0_dp, 0.0_dp, il, &
iu, -1.0_dp, m, w, cz, num_bands, cwork, rwork, iwork, ifail, info)
!~ write(stdout,*) 'w:'
!~ do n=1,ndimwin(nkp)
!~ write(stdout,'(f14.10)') w(n)
!~ enddo
!~ write(stdout,*) 'cz:'
!~ do n=1,ndimwin(nkp)
!~ write(stdout,'(6f12.8)') cz(n,il), cz(n,iu)
!~ enddo
! DEBUG
if (info .lt. 0) then
if (on_root) write (stdout, *) ' *** ERROR *** ZHPEVX WHILE DIAGONALIZING CQPQ MATRIX'
if (on_root) write (stdout, *) ' THE ', -info, ' ARGUMENT OF ZHPEVX HAD AN ILLEGAL VALUE'
call io_error('dis_proj_frozen: error')
elseif (info .gt. 0) then
if (on_root) write (stdout, *) ' *** ERROR *** ZHPEVX WHILE DIAGONALIZING CQPQ MATRIX'
if (on_root) write (stdout, *) info, 'EIGENVECTORS FAILED TO CONVERGE'
call io_error('dis_proj_frozen: error')
endif
! ENDDEBUG
! DEBUG
if (m .ne. ndimwin(nkp)) then
if (on_root) write (stdout, *) ' *** ERROR *** in dis_proj_froz'
if (on_root) write (stdout, *) ' Number of eigenvalues/vectors obtained is', &
m, ' not equal to the number asked,', ndimwin(nkp)
call io_error('dis_proj_frozen: error')
endif
! ENDDEBUG
! DEBUG
! check that the eigenvalues are between 0 and 1
if (iprint > 2) then
if (on_root) write (stdout, '(/a,i3,a,i3,a,i3,a)') ' K-point ', nkp, ' ndimwin: ', &
ndimwin(nkp), ' we want the', num_wann - ndimfroz(nkp), &
' leading eigenvector(s) of QPQ'
endif
do j = 1, ndimwin(nkp)
if (iprint > 2 .and. on_root) write (stdout, '(a,i3,a,f16.12)') ' lambda(', j, ')=', w(j)
!~[aam] if ( (w(j).lt.eps8).or.(w(j).gt.1.0_dp + eps8) ) then
if ((w(j) .lt. -eps8) .or. (w(j) .gt. 1.0_dp + eps8)) then
call io_error('dis_proj_frozen: error - Eigenvalues not between 0 and 1')
endif
enddo
! ENDDEBUG
! [ aam: sometimes the leading eigenvalues form a degenerate set that is
! of higher dimensionality than (num_wann - ndimfroz). May need to
! fix this at some point. ]
! For certain cases we have found that one of the required eigenvectors of cqpq
! has a zero eigenvalue (ie it forms a degenerate set with the frozen states).
! It depends on floating point math as whether this eigenvalue is positive
! or negative (ie +/- 1e-17). If it's positive everything is ok. If negative we
! can end up putting one of the frozen states into U_opt (and failing the later
! orthogonality check).
! This fix detects this situation. If applies we choose the eigenvectors by
! checking their orthogonality to the frozen states.
! === For version 1.0.1 we make this the default ===
if (index(devel_flag, 'no-orth-fix') == 0) then
nzero = 0; goods = 0
do j = ndimwin(nkp), ndimwin(nkp) - (num_wann - ndimfroz(nkp)) + 1, -1
if (w(j) < eps8) then
nzero = nzero + 1
else
goods = goods + 1
end if
end do
if (nzero > 0) then
if (iprint > 2 .and. on_root) then
write (stdout, *) ' '
write (stdout, '(1x,a,i0,a)') 'An eigenvalue of QPQ is close to zero at kpoint ' &
, nkp, '. Using safety check.'
write (stdout, '(1x,a,i4,a,i4)') 'We must find ', nzero, &
' eigenvectors with zero eigenvalues out of a set of ', ndimwin(nkp) - goods
endif
!First lets put the 'good' states into vamp
vmap = 0
counter = 1
do j = ndimwin(nkp), ndimwin(nkp) - goods + 1, -1
vmap(counter) = j
counter = counter + 1
end do
if (iprint > 2 .and. on_root) then
do loop_f = 1, ndimwin(nkp)
write (stdout, '(1x,a,i4,a,es13.6)') 'Eigenvector number', loop_f, ' Eigenvalue: ', w(loop_f)
do loop_v = 1, ndimwin(nkp)
write (stdout, '(20x,2f12.8)') cz(loop_v, loop_f)
end do
write (stdout, *)
end do
end if
! We need to find nzero vectors out of the remining ndimwin(nkp)-goods vectors
do loop_f = 1, nzero
do loop_v = ndimwin(nkp), 1, -1 !loop backwards for efficiency only
if (any(vmap == loop_v)) cycle
!check to see if vector is orthogonal to frozen states in u_matrix_opt
take = .true.
do m = 1, ndimfroz(nkp)
ctmp = cmplx_0
do j = 1, ndimwin(nkp)
ctmp = ctmp + conjg(u_matrix_opt(j, m, nkp))*cz(j, loop_v)
enddo
if (abs(ctmp) .gt. eps8) then
take = .false.
endif
enddo
if (take) then !vector is good and we add it to vmap
vmap(goods + loop_f) = loop_v
exit
end if
end do
end do
if (iprint > 2 .and. on_root) then
write (rep, '(i4)') num_wann - ndimfroz(nkp)
write (stdout, '(1x,a,'//trim(rep)//'(i0,1x))') 'We use the following eigenvectors: ' &
, vmap(1:(num_wann - ndimfroz(nkp)))
end if
do l = 1, num_wann - ndimfroz(nkp)
if (vmap(l) == 0) call io_error('dis_proj_froz: Ortho-fix failed to find enough vectors')
end do
! put the correct eigenvectors into u_matrix_opt, and we're all done!
counter = 1
do l = ndimfroz(nkp) + 1, num_wann
u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), vmap(counter))
counter = counter + 1
enddo
else ! we don't need to use the fix
do l = ndimfroz(nkp) + 1, num_wann
u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), il)
il = il + 1
enddo
if (il - 1 .ne. iu) then
call io_error('dis_proj_frozen: error - il-1.ne.iu (in ortho-fix)')
endif
end if
else ! if .not. using ortho-fix
! PICK THE num_wann-nDIMFROZ(NKP) LEADING EIGENVECTORS AS TRIAL STATES
! and PUT THEM RIGHT AFTER THE FROZEN STATES IN u_matrix_opt
do l = ndimfroz(nkp) + 1, num_wann
if (on_root) write (stdout, *) 'il=', il
u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), il)
il = il + 1
enddo
! DEBUG
if (il - 1 .ne. iu) then
call io_error('dis_proj_frozen: error - il-1.ne.iu')
endif
! ENDDEBUG
end if
endif ! num_wann>nDIMFROZ(NKP)
! Put the frozen states in the lowest columns of u_matrix_opt
if (ndimfroz(nkp) .gt. 0) then
do l = 1, ndimfroz(nkp)
u_matrix_opt(:, l, nkp) = cmplx_0
u_matrix_opt(indxfroz(l, nkp), l, nkp) = cmplx_1
enddo
endif
!~ write(stdout,*) 'u_matrix_opt:'
!~ do m=1,ndimwin(nkp)
!~ write(stdout,'(6f12.8)') u_matrix_opt(m,1,nkp), &
!~ u_matrix_opt(m,ndimfroz(nkp),nkp), u_matrix_opt(m,num_wann,nkp)
!~ enddo
enddo ! NKP
deallocate (cqpq, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cqpq in dis_proj_froz')
deallocate (cpq, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cpq in dis_proj_froz')
deallocate (cq_froz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cq_froz in dis_proj_froz')
deallocate (cp_s, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cp_s in dis_proj_froz')
deallocate (cz, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cz in dis_proj_froz')
deallocate (cwork, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cwork in dis_proj_froz')
deallocate (cap, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating cap in dis_proj_froz')
deallocate (rwork, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating rwork in dis_proj_froz')
deallocate (w, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating w in dis_proj_froz')
deallocate (ifail, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating ifail in dis_proj_froz')
deallocate (iwork, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating iwork in dis_proj_froz')
if (on_root) write (stdout, '(a)') ' done'
if (timing_level > 1) call io_stopwatch('dis: proj_froz', 2)
return
end subroutine dis_proj_froz