hamiltonian_get_hr Subroutine

  subroutine hamiltonian_get_hr(atom_data, dis_manifold, ham_logical, real_space_ham, &
                                print_output, ham_k, ham_r, u_matrix, u_matrix_opt, eigval, &
                                kpt_latt, real_lattice, wannier_centres, &
                                wannier_centres_translated, irvec, shift_vec, nrpts, num_bands, &
                                num_kpts, num_wann, have_disentangled, stdout, timer, error, &
                                lsitesymmetry, comm)
    !================================================!
    !
    !!  Calculate the Hamiltonian in the WF basis
    !
    !================================================!

    use w90_constants, only: cmplx_0, cmplx_i, twopi
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: atom_data_type, dis_manifold_type, print_output_type, timer_list_type
    use w90_wannier90_types, only: real_space_ham_type, ham_logical_type

    implicit none

    ! arguments
    type(ham_logical_type), intent(inout)    :: ham_logical
    type(atom_data_type), intent(in)         :: atom_data
    type(real_space_ham_type), intent(inout) :: real_space_ham
    type(print_output_type), intent(in)      :: print_output
    type(dis_manifold_type), intent(in)      :: dis_manifold
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in)           :: comm
    type(timer_list_type), intent(inout)     :: timer

    integer, intent(inout), allocatable :: shift_vec(:, :)
    integer, intent(inout)              :: irvec(:, :)
    integer, intent(inout)              :: nrpts
    integer, intent(in)                 :: num_bands
    integer, intent(in)                 :: num_kpts
    integer, intent(in)                 :: num_wann
    integer, intent(in)                 :: stdout

    real(kind=dp), intent(inout) :: wannier_centres_translated(:, :)
    real(kind=dp), intent(in)    :: real_lattice(3, 3)
    real(kind=dp), intent(in)    :: wannier_centres(:, :)
    real(kind=dp), intent(in)    :: kpt_latt(:, :)
    real(kind=dp), intent(in)    :: eigval(:, :)

    complex(kind=dp), intent(inout)              :: ham_r(:, :, :)
    complex(kind=dp), intent(in)                 :: u_matrix(:, :, :)
    complex(kind=dp), intent(in)                 :: u_matrix_opt(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: ham_k(:, :, :)

    logical, intent(in) :: lsitesymmetry  !YN:
    logical, intent(in) :: have_disentangled

    ! local variables
    integer          :: loop_kpt, i, j, m, irpt, ierr, counter
    real(kind=dp)    :: rdotk
    real(kind=dp), allocatable    :: eigval_opt(:, :) !(num_bands, num_kpts)
    real(kind=dp), allocatable    :: eigval2(:, :)    !(num_wann, num_kpts)
    real(kind=dp)    :: irvec_tmp(3)
    complex(kind=dp), allocatable :: utmp(:, :)       !(num_bands, num_wann)
    complex(kind=dp) :: fac

    if (print_output%timing_level > 1) call io_stopwatch_start('hamiltonian: get_hr', timer)

    if (ham_logical%have_ham_r) then
      if (ham_logical%have_translated .eqv. ham_logical%use_translation) then
        goto 200
      else
        goto 100
      end if
    end if

    if (ham_logical%have_ham_k) go to 100

    ham_k = cmplx_0

    allocate (eigval2(num_wann, num_kpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating eigval2 in hamiltonian_get_hr', comm)
      return
    end if

    eigval2 = 0.0_dp

    if (have_disentangled) then

      ! start allocation of eigval_opt, utmp; used only if have_disentangled.
      allocate (eigval_opt(num_bands, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating eigval_opt in hamiltonian_get_hr', comm)
        return
      end if

      allocate (utmp(num_bands, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating utmp in hamiltonian_get_hr', comm)
        return
      end if

      eigval_opt = 0.0_dp
      ! end allocation of eigval_opt, utmp

      ! slim down eigval to contain states within the outer window

      do loop_kpt = 1, num_kpts
        counter = 0
        do j = 1, num_bands
          if (dis_manifold%lwindow(j, loop_kpt)) then
            counter = counter + 1
            eigval_opt(counter, loop_kpt) = eigval(j, loop_kpt)
          end if
        end do
      end do

      ! rotate eigval into the optimal subspace
      ! in general eigval would be a matrix at each kpoints
      ! but we choose u_matrix_opt such that the Hamiltonian is
      ! diagonal at each kpoint. (I guess we should check it here)

      if (.not. lsitesymmetry) then
        do loop_kpt = 1, num_kpts
          do j = 1, num_wann
            do m = 1, dis_manifold%ndimwin(loop_kpt)
              eigval2(j, loop_kpt) = eigval2(j, loop_kpt) + eigval_opt(m, loop_kpt)* &
                                     real(conjg(u_matrix_opt(m, j, loop_kpt))* &
                                          u_matrix_opt(m, j, loop_kpt), dp)
            end do
          end do
        end do
      else
        ! u_matrix_opt are not the eigenvectors of the Hamiltonian any more
        ! so we have to calculate ham_k in the following way
        do loop_kpt = 1, num_kpts
          utmp(1:dis_manifold%ndimwin(loop_kpt), :) = &
            matmul(u_matrix_opt(1:dis_manifold%ndimwin(loop_kpt), :, loop_kpt), &
                   u_matrix(:, :, loop_kpt))
          do j = 1, num_wann
            do i = 1, j
              do m = 1, dis_manifold%ndimwin(loop_kpt)
                ham_k(i, j, loop_kpt) = ham_k(i, j, loop_kpt) + eigval_opt(m, loop_kpt)* &
                                        conjg(utmp(m, i))*utmp(m, j)
              end do
              if (i .lt. j) ham_k(j, i, loop_kpt) = conjg(ham_k(i, j, loop_kpt))
            end do
          end do
        end do
      end if

    else
      eigval2(1:num_wann, :) = eigval(1:num_wann, :)
    end if

    ! At this point eigval2 contains num_wann values which belong to the wannier subspace.

    ! Rotate Hamiltonian into the basis of smooth bloch states
    !          H(k)=U^{dagger}(k).H_0(k).U(k)
    ! Note: we enforce hermiticity here

    if (.not. lsitesymmetry .or. .not. have_disentangled) then
      do loop_kpt = 1, num_kpts
        do j = 1, num_wann
          do i = 1, j
            do m = 1, num_wann
              ham_k(i, j, loop_kpt) = ham_k(i, j, loop_kpt) + eigval2(m, loop_kpt)* &
                                      conjg(u_matrix(m, i, loop_kpt))*u_matrix(m, j, loop_kpt)
            end do
            if (i .lt. j) ham_k(j, i, loop_kpt) = conjg(ham_k(i, j, loop_kpt))
          end do
        end do
      end do
    end if

    ham_logical%have_ham_k = .true.

100 continue

    ! Fourier transform rotated hamiltonian into WF basis
    ! H_ij(k) --> H_ij(R) = (1/N_kpts) sum_k e^{-ikR} H_ij(k)
!~    if (.not.allocated(ham_r)) then
!~      allocate(ham_r(num_wann,num_wann,nrpts),stat=ierr)
!~      if (ierr/=0) call io_error('Error in allocating ham_r in hamiltonian_get_hr')
!~    end if

    ham_r = cmplx_0

    if (.not. ham_logical%use_translation) then

      do irpt = 1, nrpts
        do loop_kpt = 1, num_kpts
          rdotk = twopi*dot_product(kpt_latt(:, loop_kpt), real(irvec(:, irpt), dp))
          fac = exp(-cmplx_i*rdotk)/real(num_kpts, dp)
          ham_r(:, :, irpt) = ham_r(:, :, irpt) + fac*ham_k(:, :, loop_kpt)
        end do
      end do

      ham_logical%have_translated = .false.

    else

      allocate (shift_vec(3, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating shift_vec in hamiltonian_get_hr', comm)
        return
      end if
      call internal_translate_centres(atom_data, real_space_ham, real_lattice, wannier_centres, &
                                      wannier_centres_translated, shift_vec, print_output%iprint, &
                                      num_wann, error)
      if (allocated(error)) return

      do irpt = 1, nrpts
        do loop_kpt = 1, num_kpts
          do i = 1, num_wann
            do j = 1, num_wann
              ! ham_r(j,i,irpt)
              ! interaction btw j at 0 and i at irvec(:,irpt)
              irvec_tmp(:) = irvec(:, irpt) + shift_vec(:, i) - shift_vec(:, j)
              rdotk = twopi*dot_product(kpt_latt(:, loop_kpt), real(irvec_tmp(:), dp))
              fac = exp(-cmplx_i*rdotk)/real(num_kpts, dp)
              ham_r(j, i, irpt) = ham_r(j, i, irpt) + fac*ham_k(j, i, loop_kpt)
            end do
          end do
        end do
      end do

      ham_logical%have_translated = .true.

    end if

    ! [lp] if required, compute the minimum diistances
!     if (use_ws_distance) then
!         allocate(irdist_ws(3,ndegenx,num_wann,num_wann,nrpts),stat=ierr)
!         if (ierr/=0) call io_error('Error in allocating irdist_ws in hamiltonian_get_hr')
!         allocate(wdist_ndeg(num_wann,num_wann,nrpts),stat=ierr)
!         if (ierr/=0) call io_error('Error in allocating wcenter_ndeg in hamiltonian_get_hr')
    !
!         call ws_translate_dist(nrpts, irvec)
!     endif

    ham_logical%have_ham_r = .true.

200 continue

    if (allocated(shift_vec)) then
      deallocate (shift_vec, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating shift_vec in hamiltonian_get_hr', comm)
        return
      end if
    end if

    if (allocated(eigval2)) then
      deallocate (eigval2, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating eigval2 in hamiltonian_get_hr', comm)
        return
      end if
    end if

    if (allocated(eigval_opt)) then
      deallocate (eigval_opt, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating eigval_opt in hamiltonian_get_hr', comm)
        return
      end if
    end if

    if (allocated(utmp)) then
      deallocate (utmp, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating utmp in hamiltonian_get_hr', comm)
        return
      end if
    end if

    if (print_output%timing_level > 1) call io_stopwatch_stop('hamiltonian: get_hr', timer)

    return

  contains

    !================================================!
    subroutine internal_translate_centres(atom_data, real_space_ham, real_lattice, &
                                          wannier_centres, wannier_centres_translated, shift_vec, &
                                          iprint, num_wann, error)
      !================================================!
      !
      !! Translate the centres of the WF into the home cell
      !
      !================================================!

      use w90_utility, only: utility_cart_to_frac, utility_frac_to_cart, utility_inverse_mat
      use w90_types, only: atom_data_type
      use w90_wannier90_types, only: real_space_ham_type

      implicit none

      ! arguments
      type(atom_data_type), intent(in) :: atom_data
      type(real_space_ham_type), intent(inout) :: real_space_ham
      type(w90_error_type), allocatable, intent(out) :: error

      integer, intent(inout) :: shift_vec(:, :)
      integer, intent(in)    :: iprint
      integer, intent(in)    :: num_wann

      real(kind=dp), intent(inout) :: wannier_centres_translated(:, :)
      real(kind=dp), intent(in)    :: real_lattice(3, 3)
      real(kind=dp), intent(in)    :: wannier_centres(:, :)

      ! local variables
      integer :: iw, ierr, nat, nsp, ind
      real(kind=dp)              :: inv_lattice(3, 3)
      real(kind=dp), allocatable :: r_home(:, :), r_frac(:, :)
      real(kind=dp)              :: c_pos_cart(3), c_pos_frac(3)
      real(kind=dp)              :: r_frac_min(3)

!~      if (.not.allocated(wannier_centres_translated)) then
!~         allocate(wannier_centres_translated(3,num_wann),stat=ierr)
!~         if (ierr/=0) call io_error('Error in allocating wannier_centres_translated &
!~              &in internal_translate_wannier_centres')
!~      end if

      allocate (r_home(3, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating r_home in internal_translate_centres', comm)
        return
      end if
      allocate (r_frac(3, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating r_frac in internal_translate_centres', comm)
        return
      end if
      r_home = 0.0_dp; r_frac = 0.0_dp

      call utility_inverse_mat(real_lattice, inv_lattice)
      if (real_space_ham%automatic_translation) then
        ! Calculate centre of atomic positions
        c_pos_cart = 0.0_dp; c_pos_frac = 0.0_dp
        do nsp = 1, atom_data%num_species
          do nat = 1, atom_data%species_num(nsp)
            c_pos_cart(:) = c_pos_cart(:) + atom_data%pos_cart(:, nat, nsp)
          end do
        end do
        c_pos_cart = c_pos_cart/atom_data%num_atoms
        ! Cartesian --> fractional
        call utility_cart_to_frac(c_pos_cart, real_space_ham%translation_centre_frac, inv_lattice)
      end if
      ! Wannier function centres will be in [c_pos_frac-0.5,c_pos_frac+0.5]
      r_frac_min(:) = real_space_ham%translation_centre_frac(:) - 0.5_dp

      ! Cartesian --> fractional
      do iw = 1, num_wann
        call utility_cart_to_frac(wannier_centres(:, iw), r_frac(:, iw), inv_lattice)
        ! Rationalise r_frac - r_frac_min to interval [0,1]
        !  by applying shift of -floor(r_frac - r_frac_min)
        shift_vec(:, iw) = -floor(r_frac(:, iw) - r_frac_min(:))
        r_frac(:, iw) = r_frac(:, iw) + real(shift_vec(:, iw), dp)
        ! Fractional --> Cartesian
        call utility_frac_to_cart(r_frac(:, iw), r_home(:, iw), real_lattice)
      end do

      ! NEVER overwrite wannier_centres
      !wannier_centres = r_home

      if (iprint > 0) then
        write (stdout, '(1x,a)') 'Translated centres'
        write (stdout, '(4x,a,3f10.6)') 'translation centre in fractional coordinate:', &
          real_space_ham%translation_centre_frac(:)
        do iw = 1, num_wann
          write (stdout, 888) iw, (r_home(ind, iw)*print_output%lenconfac, ind=1, 3)
        end do
        write (stdout, '(1x,a78)') repeat('-', 78)
        write (stdout, *)
      end if
      wannier_centres_translated = r_home

      deallocate (r_frac, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating r_frac in internal_translate_centres', comm)
        return
      end if
      deallocate (r_home, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating r_home in internal_translate_centres', comm)
        return
      end if

      return

888   format(2x, 'WF centre ', i5, 2x, '(', f10.6, ',', f10.6, ',', f10.6, ' )')

    end subroutine internal_translate_centres

  end subroutine hamiltonian_get_hr