get_HH_R Subroutine

  subroutine get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, &
                      v_matrix, eigval, real_lattice, scissors_shift, num_bands, num_kpts, &
                      num_wann, num_valence_bands, effective_model, have_disentangled, seedname, &
                      ws_distance, ws_region, stdout, timer, error, comm)
    !================================================
    !
    !! computes <0n|H|Rm>, in eV
    !! (pwscf uses Ry, but pw2wannier90 converts to eV)
    !
    !================================================

    use w90_postw90_types, only: wigner_seitz_type
    use w90_types, only: dis_manifold_type, print_output_type, timer_list_type, &
                         ws_distance_type, ws_region_type
    implicit none

    ! arguments
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(print_output_type), intent(in) :: print_output
    type(w90_comm_type), intent(in) :: comm
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(in) :: ws_distance
    type(ws_region_type), intent(in) :: ws_region
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

    integer, intent(in) :: num_bands, num_kpts, num_wann, num_valence_bands, stdout

    real(kind=dp), intent(in) :: eigval(:, :), real_lattice(3, 3)
    real(kind=dp), intent(in) :: scissors_shift
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable :: HH_R_temp(:, :, :)

    character(len=50), intent(in) :: seedname

    logical, intent(in) :: have_disentangled
    logical, intent(in) :: effective_model

    ! local variables
    integer                       :: i, j, n, m, ii, ik, winmin_q, file_unit, &
                                     ir, io, idum, ivdum(3), ivdum_old(3)
    integer, allocatable          :: num_states(:)
    real(kind=dp)                 :: rdum_real, rdum_imag
    complex(kind=dp), allocatable :: HH_q(:, :, :)
    logical                       :: new_ir

    complex(kind=dp), allocatable :: sciss_q(:, :, :)
    complex(kind=dp), allocatable :: sciss_R(:, :, :)

    logical :: on_root = .false.
    if (mpirank(comm) == 0) on_root = .true.

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_start('get_oper: get_HH_R', timer)

    if (wigner_seitz%nrpts < 1) then
      call set_error_fatal(error, 'Error: wigner_setiz%nrpts incorrect at get_HH_R() ', comm)
      return
    end if

    allocate (HH_R_temp(num_wann, num_wann, wigner_seitz%nrpts))

    if (.not. allocated(HH_R)) then
      allocate (HH_R(num_wann, num_wann, wigner_seitz%nrpts_pw90))
    else
      if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
        call io_stopwatch_stop('get_oper: get_HH_R', timer)
      return
    end if

    ! Real-space Hamiltonian H(R) is read from file
    !
    if (effective_model) then
      HH_R = cmplx_0
      if (on_root) then
        write (stdout, '(/a)') ' Reading real-space Hamiltonian from file ' &
          //trim(seedname)//'_HH_R.dat'
        open (newunit=file_unit, file=trim(seedname)//'_HH_R.dat', form='formatted', &
              status='old', err=101)
        read (file_unit, *) ! header
        read (file_unit, *) idum ! num_wann
        read (file_unit, *) idum ! nrpts
        ir = 1
        new_ir = .true.
        ivdum_old(:) = 0
        n = 1
        do
          read (file_unit, '(5I5,2F12.6)', iostat=io) ivdum(1:3), j, i, &
            rdum_real, rdum_imag
          if (io < 0) exit ! reached end of file
          if (i < 1 .or. i > num_wann .or. j < 1 .or. j > num_wann) then
            write (stdout, *) 'num_wann=', num_wann, '  i=', i, '  j=', j
            call set_error_fatal(error, 'Error in get_HH_R: orbital indices out of bounds', comm)
            return
          end if
          if (n > 1) then
            if (ivdum(1) /= ivdum_old(1) .or. ivdum(2) /= ivdum_old(2) .or. &
                ivdum(3) /= ivdum_old(3)) then
              ir = ir + 1
              new_ir = .true.
            else
              new_ir = .false.
            end if
          end if
          ivdum_old = ivdum
          ! Note that the same (j,i,ir) may occur more than once in
          ! the file seedname_HH_R.dat, hence the addition instead
          ! of a simple equality. (This has to do with the way the
          ! Berlijn effective Hamiltonian algorithm is
          ! implemented.)
          HH_R(j, i, ir) = HH_R(j, i, ir) + cmplx(rdum_real, rdum_imag, kind=dp)
          if (new_ir) then
            wigner_seitz%irvec(:, ir) = ivdum(:)
            if (ivdum(1) == 0 .and. ivdum(2) == 0 .and. ivdum(3) == 0) wigner_seitz%rpt_origin = ir
          end if
          n = n + 1
        end do
        close (file_unit)
        if (ir /= wigner_seitz%nrpts) then
          write (stdout, *) 'ir=', ir, '  nrpts=', wigner_seitz%nrpts
          call set_error_fatal(error, 'Error in get_HH_R: inconsistent nrpts values', comm)
          return
        end if
        do ir = 1, wigner_seitz%nrpts
          wigner_seitz%crvec(:, ir) = matmul(transpose(real_lattice), wigner_seitz%irvec(:, ir))
        end do
        wigner_seitz%ndegen(:) = 1 ! This is assumed when reading HH_R from file
        !
        wigner_seitz%nrpts_pw90 = wigner_seitz%nrpts
        wigner_seitz%irvec_pw90 = wigner_seitz%irvec
        wigner_seitz%crvec_pw90 = wigner_seitz%crvec
        !
        ! TODO: Implement scissors in this case? Need to choose a
        ! uniform k-mesh (the scissors correction is applied in
        ! k-space) and then proceed as below, Fourier transforming
        ! back to real space and adding to HH_R_temp, Hopefully the
        ! result converges (rapidly) with the k-mesh density, but
        ! one should check
        !
        if (abs(scissors_shift) > 1.0e-7_dp) then
          call set_error_input(error, 'Error in get_HH_R: scissors shift not implemented for ' &
                               //'effective_model=T', comm)
          return
        end if
      end if
      call comms_bcast(HH_R(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, error, comm)
      if (allocated(error)) return
      call comms_bcast(wigner_seitz%ndegen(1), wigner_seitz%nrpts, error, comm)
      if (allocated(error)) return
      call comms_bcast(wigner_seitz%irvec(1, 1), 3*wigner_seitz%nrpts, error, comm)
      if (allocated(error)) return
      call comms_bcast(wigner_seitz%crvec(1, 1), 3*wigner_seitz%nrpts, error, comm)
      if (allocated(error)) return
      if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
        call io_stopwatch_stop('get_oper: get_HH_R', timer)
      return
    end if

    ! Everything below is only executed if effective_model==False (default)

    ! Real-space Hamiltonian H(R) is calculated by Fourier
    ! transforming H(q) defined on the ab-initio reciprocal mesh
    !
    allocate (HH_q(num_wann, num_wann, num_kpts))
    allocate (num_states(num_kpts))

    HH_q = cmplx_0
    do ik = 1, num_kpts
      if (have_disentangled) then
        num_states(ik) = dis_manifold%ndimwin(ik)
      else
        num_states(ik) = num_wann
      end if

      call get_win_min(num_bands, dis_manifold, ik, winmin_q, have_disentangled)
      do m = 1, num_wann
        do n = 1, m
          do i = 1, num_states(ik)
            ii = winmin_q + i - 1
            HH_q(n, m, ik) = HH_q(n, m, ik) &
                             + conjg(v_matrix(i, n, ik))*eigval(ii, ik) &
                             *v_matrix(i, m, ik)
          end do
          HH_q(m, n, ik) = conjg(HH_q(n, m, ik))
        end do
      end do
    end do

    call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, HH_q, HH_R_temp)

    ! Scissors correction for an insulator: shift conduction bands upwards by
    ! scissors_shift eV
    !
    if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then
      allocate (sciss_R(num_wann, num_wann, wigner_seitz%nrpts))
      allocate (sciss_q(num_wann, num_wann, num_kpts))
      sciss_q = cmplx_0
      do ik = 1, num_kpts
        do j = 1, num_wann
          do i = 1, j
            do m = 1, num_valence_bands
              sciss_q(i, j, ik) = sciss_q(i, j, ik) - &
                                  conjg(u_matrix(m, i, ik))*u_matrix(m, j, ik)
            end do
            sciss_q(j, i, ik) = conjg(sciss_q(i, j, ik))
          end do
        end do
      end do

      call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, sciss_q, &
                          sciss_R)
      do n = 1, num_wann
        sciss_R(n, n, wigner_seitz%rpt_origin) = sciss_R(n, n, wigner_seitz%rpt_origin) + 1.0_dp
      end do
      sciss_R = sciss_R*scissors_shift
      HH_R_temp = HH_R_temp + sciss_R
    end if

    ! Apply degeneracy factor and reorder according to the wigner-seitz vectors
    call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, HH_R_temp, HH_R)

    if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
      call io_stopwatch_stop('get_oper: get_HH_R', timer)
    return

101 call set_error_file(error, 'Error in get_HH_R: problem opening file '// &
                        trim(seedname)//'_HH_R.dat', comm)
    return

  end subroutine get_HH_R