get_HH_R Subroutine

public subroutine get_HH_R()

Uses

  • proc~~get_hh_r~~UsesGraph proc~get_hh_r get_HH_R module~w90_constants w90_constants proc~get_hh_r->module~w90_constants module~w90_io w90_io proc~get_hh_r->module~w90_io module~w90_postw90_common w90_postw90_common proc~get_hh_r->module~w90_postw90_common module~w90_parameters w90_parameters proc~get_hh_r->module~w90_parameters module~w90_comms w90_comms proc~get_hh_r->module~w90_comms module~w90_io->module~w90_constants module~w90_postw90_common->module~w90_constants module~w90_postw90_common->module~w90_comms module~w90_parameters->module~w90_constants module~w90_parameters->module~w90_io module~w90_comms->module~w90_constants module~w90_comms->module~w90_io

computes <0n|H|Rm>, in eV (pwscf uses Ry, but pw2wannier90 converts to eV)

Arguments

None

Calls

proc~~get_hh_r~~CallsGraph proc~get_hh_r get_HH_R proc~io_error io_error proc~get_hh_r->proc~io_error proc~get_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~io_file_unit io_file_unit proc~get_hh_r->proc~io_file_unit proc~fourier_q_to_r fourier_q_to_R proc~get_hh_r->proc~fourier_q_to_r

Called by

proc~~get_hh_r~~CalledByGraph proc~get_hh_r get_HH_R proc~gyrotropic_main gyrotropic_main proc~gyrotropic_main->proc~get_hh_r proc~gyrotropic_get_k_list gyrotropic_get_k_list proc~gyrotropic_main->proc~gyrotropic_get_k_list proc~k_path k_path proc~k_path->proc~get_hh_r proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~k_path->proc~berry_get_imfgh_klist proc~berry_get_imf_klist berry_get_imf_klist proc~k_path->proc~berry_get_imf_klist proc~wham_get_eig_uu_hh_aa_sc_tb_conv wham_get_eig_UU_HH_AA_sc_TB_conv proc~wham_get_eig_uu_hh_aa_sc_tb_conv->proc~get_hh_r proc~wham_get_eig_deleig wham_get_eig_deleig proc~wham_get_eig_deleig->proc~get_hh_r proc~berry_main berry_main proc~berry_main->proc~get_hh_r proc~berry_get_kubo_k berry_get_kubo_k proc~berry_main->proc~berry_get_kubo_k proc~berry_main->proc~berry_get_imfgh_klist proc~berry_get_sc_klist berry_get_sc_klist proc~berry_main->proc~berry_get_sc_klist proc~berry_main->proc~berry_get_imf_klist proc~k_slice k_slice proc~k_slice->proc~get_hh_r proc~k_slice->proc~wham_get_eig_deleig proc~k_slice->proc~berry_get_imfgh_klist proc~k_slice->proc~berry_get_imf_klist proc~wham_get_eig_uu_hh_jjlist wham_get_eig_UU_HH_JJlist proc~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~geninterp_main geninterp_main proc~geninterp_main->proc~get_hh_r proc~spin_get_moment spin_get_moment proc~spin_get_moment->proc~get_hh_r proc~dos_main dos_main proc~dos_main->proc~get_hh_r proc~dos_main->proc~wham_get_eig_deleig proc~wham_get_eig_uu_hh_aa_sc wham_get_eig_UU_HH_AA_sc proc~wham_get_eig_uu_hh_aa_sc->proc~get_hh_r proc~calctdfanddos calcTDFandDOS proc~calctdfanddos->proc~get_hh_r proc~calctdfanddos->proc~wham_get_eig_deleig proc~berry_get_kubo_k->proc~wham_get_eig_deleig proc~berry_get_imfgh_klist->proc~wham_get_eig_uu_hh_jjlist proc~gyrotropic_get_k_list->proc~wham_get_eig_deleig proc~gyrotropic_get_k_list->proc~berry_get_imfgh_klist proc~gyrotropic_get_k_list->proc~berry_get_imf_klist proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc_tb_conv proc~berry_get_sc_klist->proc~wham_get_eig_deleig proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc proc~boltzwann_main boltzwann_main proc~boltzwann_main->proc~calctdfanddos proc~berry_get_imf_klist->proc~berry_get_imfgh_klist

Contents

Source Code


Source Code

  subroutine get_HH_R
    !======================================================
    !
    !! computes <0n|H|Rm>, in eV
    !! (pwscf uses Ry, but pw2wannier90 converts to eV)
    !
    !======================================================

    use w90_constants, only: dp, cmplx_0
    use w90_io, only: io_error, stdout, io_stopwatch, &
      io_file_unit, seedname
    use w90_parameters, only: num_wann, ndimwin, num_kpts, num_bands, &
      eigval, u_matrix, have_disentangled, &
      timing_level, scissors_shift, &
      num_valence_bands, effective_model, &
      real_lattice
    use w90_postw90_common, only: nrpts, rpt_origin, v_matrix, ndegen, irvec, crvec
    use w90_comms, only: on_root, comms_bcast

    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

    !ivo
    complex(kind=dp), allocatable :: sciss_q(:, :, :)
    complex(kind=dp), allocatable :: sciss_R(:, :, :)
    real(kind=dp)                 :: sciss_shift

    if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 1)

    if (.not. allocated(HH_R)) then
      allocate (HH_R(num_wann, num_wann, nrpts))
    else
      if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
      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'
        file_unit = io_file_unit()
        open (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 io_error &
              ('Error in get_HH_R: orbital indices out of bounds')
          endif
          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.
            endif
          endif
          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
            irvec(:, ir) = ivdum(:)
            if (ivdum(1) == 0 .and. ivdum(2) == 0 .and. ivdum(3) == 0) rpt_origin = ir
          endif
          n = n + 1
        enddo
        close (file_unit)
        if (ir /= nrpts) then
          write (stdout, *) 'ir=', ir, '  nrpts=', nrpts
          call io_error('Error in get_HH_R: inconsistent nrpts values')
        endif
        do ir = 1, nrpts
          crvec(:, ir) = matmul(transpose(real_lattice), irvec(:, ir))
        end do
        ndegen(:) = 1 ! This is assumed when reading HH_R from file
        !
        ! 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, Hopefully the
        ! result converges (rapidly) with the k-mesh density, but
        ! one should check
        !
        if (abs(scissors_shift) > 1.0e-7_dp) &
          call io_error( &
          'Error in get_HH_R: scissors shift not implemented for ' &
          //'effective_model=T')
      endif
      call comms_bcast(HH_R(1, 1, 1), num_wann*num_wann*nrpts)
      call comms_bcast(ndegen(1), nrpts)
      call comms_bcast(irvec(1, 1), 3*nrpts)
      call comms_bcast(crvec(1, 1), 3*nrpts)
      if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
      return
    endif

    ! 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) = ndimwin(ik)
      else
        num_states(ik) = num_wann
      endif
      call get_win_min(ik, winmin_q)
      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)
          enddo
          HH_q(m, n, ik) = conjg(HH_q(n, m, ik))
        enddo
      enddo
    enddo
    call fourier_q_to_R(HH_q, HH_R)

    ! 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, 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)
            enddo
            sciss_q(j, i, ik) = conjg(sciss_q(i, j, ik))
          enddo
        enddo
      enddo
      call fourier_q_to_R(sciss_q, sciss_R)
      do n = 1, num_wann
        sciss_R(n, n, rpt_origin) = sciss_R(n, n, rpt_origin) + 1.0_dp
      end do
      sciss_R = sciss_R*scissors_shift
      HH_R = HH_R + sciss_R
    endif

    if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
    return

101 call io_error('Error in get_HH_R: problem opening file '// &
                  trim(seedname)//'_HH_R.dat')

  end subroutine get_HH_R