sitesym_dis_extract_symmetry Subroutine

public subroutine sitesym_dis_extract_symmetry(ik, n, zmat, lambda, umat)

Uses

  • proc~~sitesym_dis_extract_symmetry~~UsesGraph proc~sitesym_dis_extract_symmetry sitesym_dis_extract_symmetry module~w90_parameters w90_parameters proc~sitesym_dis_extract_symmetry->module~w90_parameters module~w90_constants w90_constants module~w90_parameters->module~w90_constants module~w90_io w90_io module~w90_parameters->module~w90_io module~w90_io->module~w90_constants

Arguments

Type IntentOptional AttributesName
integer, intent(in) :: ik
integer, intent(in) :: n
complex(kind=dp), intent(in) :: zmat(num_bands,num_bands,num_kpts)
complex(kind=dp), intent(out) :: lambda(num_wann,num_wann)
complex(kind=dp), intent(inout) :: umat(num_bands,num_wann,num_kpts)

Calls

proc~~sitesym_dis_extract_symmetry~~CallsGraph proc~sitesym_dis_extract_symmetry sitesym_dis_extract_symmetry proc~io_error io_error proc~sitesym_dis_extract_symmetry->proc~io_error zhpgvx zhpgvx proc~sitesym_dis_extract_symmetry->zhpgvx

Contents


Source Code

  subroutine sitesym_dis_extract_symmetry(ik, n, zmat, lambda, umat)
    !==================================================================!
    !                                                                  !
    !   minimize Omega_I by steepest descendent                        !
    !                                                                  !
    !   delta U_{mu I}(k) = Z_{mu mu'}*U_{mu' I}(k)                    !
    !                        - \sum_{J} lambda_{JI} U_{mu J}(k)        !
    !   lambda_{JI}=U^{*}_{mu J} Z_{mu mu'} U_{mu' I}                  !
    !                                                                  !
    !==================================================================!
    use w90_parameters, only: num_bands, num_wann, num_kpts

    implicit none

    integer, intent(in) :: ik, n
    complex(kind=dp), intent(in) :: zmat(num_bands, num_bands, num_kpts)
    complex(kind=dp), intent(out) :: lambda(num_wann, num_wann)
    complex(kind=dp), intent(inout) :: umat(num_bands, num_wann, num_kpts)

    complex(kind=dp) :: umatnew(num_bands, num_wann)
    complex(kind=dp) :: ZU(num_bands, num_wann)
    complex(kind=dp) :: deltaU(num_bands, num_wann), carr(num_bands)
    integer :: i, m, INFO, IFAIL(2), IWORK(5*2)
    complex(kind=dp) :: HP(3), SP(3), V(2, 2), CWORK(2*2)
    real(kind=dp)    :: W(2), RWORK(7*2), sp3
    integer :: iter
    integer, parameter :: niter = 50

    do iter = 1, niter
      !  Z*U
      call zgemm('N', 'N', n, num_wann, n, cmplx_1, &
                 zmat(:, :, ik), num_bands, umat(:, :, ik), num_bands, &
                 cmplx_0, ZU, num_bands)
      ! lambda = U^{+}*Z*U
      call zgemm('C', 'N', num_wann, num_wann, n, cmplx_1, &
                 umat(:, :, ik), num_bands, ZU, num_bands, &
                 cmplx_0, lambda, num_wann)

      deltaU(:, :) = ZU(:, :) - matmul(umat(:, :, ik), lambda)
      if (sum(abs(deltaU(:n, :))) .lt. 1e-10) return
      ! band-by-band minimization
      do i = 1, num_wann
        ! diagonalize 2x2 matrix
        HP(1) = real(dot_product(umat(1:n, i, ik), ZU(1:n, i)), kind=dp)
        HP(2) = dot_product(ZU(1:n, i), deltaU(1:n, i)) ! (1,2) matrix element
        carr(1:n) = matmul(zmat(1:n, 1:n, ik), deltaU(1:n, i))
        HP(3) = real(dot_product(deltaU(1:n, i), carr(1:n)), kind=dp) ! (2,2)

        SP(1) = real(dot_product(umat(1:n, i, ik), umat(1:n, i, ik)), kind=dp)
        SP(2) = dot_product(umat(1:n, i, ik), deltaU(1:n, i))
        SP(3) = real(dot_product(deltaU(1:n, i), deltaU(1:n, i)), kind=dp)

        sp3 = real(SP(3), kind=dp)
        if (abs(sp3) .lt. 1e-10) then
          umatnew(:, i) = umat(:, i, ik)
          cycle
        endif
        call ZHPGVX(1, 'V', 'A', 'U', 2, HP, SP, 0.0_dp, 0.0_dp, 0, 0, &
                    -1.0_dp, m, W, V, 2, CWORK, RWORK, IWORK, IFAIL, INFO)
        if (INFO .ne. 0) then
          write (stdout, *) 'error in sitesym_dis_extract_symmetry: INFO=', INFO
          if (INFO .gt. 0) then
            if (INFO .le. 2) then
              write (stdout, *) INFO, ' eigenvectors failed to converge'
              write (stdout, *) IFAIL(1:INFO)
            else
              write (stdout, *) ' S is not positive definite'
              write (stdout, *) 'sp3=', sp3
            endif
            call io_error('error at sitesym_dis_extract_symmetry')
          endif
        endif
        ! choose the larger eigenstate
        umatnew(:, i) = V(1, 2)*umat(:, i, ik) + V(2, 2)*deltaU(:, i)
      enddo ! i
      call symmetrize_ukirr(ik2ir(ik), num_bands, umatnew, n)
      umat(:, :, ik) = umatnew(:, :)
    enddo ! iter

    return
  end subroutine sitesym_dis_extract_symmetry