internal_zmatrix Subroutine

  subroutine internal_zmatrix(cbw, cmtrx, m_matrix_orig_local, u_matrix_opt, wb, indxnfroz, &
                              ndimfroz, ndimwin, nnlist, nkp, nkp_loc, nntot, num_bands, num_wann, &
                              timing_level, timer)
    !================================================!
    !
    !! Compute the Z-matrix
    !
    !================================================!
    use w90_constants, only: dp, cmplx_0, cmplx_1
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: timer_list_type

    implicit none

    ! arguments
    type(timer_list_type), intent(inout) :: timer

    integer, intent(in) :: num_bands, num_wann
    integer, intent(in) :: timing_level
    integer, intent(in) :: ndimwin(:)
    integer, intent(in) :: nntot, nnlist(:, :)
    integer, intent(in) :: indxnfroz(:, :)
    integer, intent(in) :: ndimfroz(:)
    integer, intent(in) :: nkp
    integer, intent(in) :: nkp_loc

    real(kind=dp), intent(in) :: wb(:)

    complex(kind=dp), intent(in) :: u_matrix_opt(:, :, :)
    complex(kind=dp), intent(in) :: m_matrix_orig_local(:, :, :, :)
    complex(kind=dp), intent(inout) :: cbw(:, :)
    complex(kind=dp), intent(out) :: cmtrx(:, :)
    !! (M,N)-TH ENTRY IN THE (NDIMWIN(NKP)-NDIMFROZ(NKP)) x (NDIMWIN(NKP)-NDIMFRO
    !! HERMITIAN MATRIX AT THE NKP-TH K-POINT

    ! local variables
    integer          :: l, m, n, p, q, nn, nkp2, ndimk
    complex(kind=dp) :: csum

    if (timing_level > 1) call io_stopwatch_start('dis: extract: zmatrix', timer)

    cmtrx = cmplx_0
    ndimk = ndimwin(nkp) - ndimfroz(nkp)
    do nn = 1, nntot
      nkp2 = nnlist(nkp, nn)
      call zgemm('N', 'N', num_bands, num_wann, ndimwin(nkp2), cmplx_1, &
                 m_matrix_orig_local(:, :, nn, nkp_loc), num_bands, u_matrix_opt(:, :, nkp2), &
                 num_bands, cmplx_0, cbw, num_bands)
      do n = 1, ndimk
        q = indxnfroz(n, nkp)
        do m = 1, n
          p = indxnfroz(m, nkp)
          csum = cmplx_0
          do l = 1, num_wann
            csum = csum + cbw(p, l)*conjg(cbw(q, l))
          end do
          cmtrx(m, n) = cmtrx(m, n) + cmplx(wb(nn), 0.0_dp, kind=dp)*csum
          cmtrx(n, m) = conjg(cmtrx(m, n))
        end do
      end do
    end do

    if (timing_level > 1) call io_stopwatch_stop('dis: extract: zmatrix', timer)

    return
    !================================================!
  end subroutine internal_zmatrix