get_FF_R Subroutine

  subroutine get_FF_R(num_bands, num_kpts, num_wann, wigner_seitz, ws_distance, ws_region, v_matrix, &
                      FF_R, dis_manifold, kmesh_info, kpt_latt, print_output, have_disentangled, stdout, &
                      seedname, timer, error, comm)
    !================================================
    !
    !! FF_ab(R) = <0|r_a.(r-R)_b|R> is the Fourier transform of
    !! FF_ab(k) = <del_a u|del_b u> (a=alpha,b=beta)
    !
    !================================================
    use w90_postw90_types, only: wigner_seitz_type
    use w90_types, only: dis_manifold_type, kmesh_info_type, ws_distance_type, ws_region_type, &
                         print_output_type, timer_list_type

    implicit none

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

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

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

    complex(kind=dp), intent(in) :: v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: FF_R(:, :, :, :, :) ! <0|r_alpha.(r-R)_beta|R>
    complex(kind=dp), allocatable :: FF_R_temp(:, :, :, :, :)

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

    logical, intent(in) :: have_disentangled

    ! local variables
    integer          :: i, j, ii, jj, m, n, a, b, nn1, nn2, ik, nb_tmp, nkp_tmp, nntot_tmp, &
                        uIu_in, qb1, qb2, winmin_qb1, winmin_qb2

    integer, allocatable          :: num_states(:)
    complex(kind=dp), allocatable :: FF_q(:, :, :, :, :)
    complex(kind=dp), allocatable :: Lo_qb1_q_qb2(:, :)
    complex(kind=dp), allocatable :: L_qb1_q_qb2(:, :)
    character(len=60)             :: header
    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_FF_R', timer)

    allocate (FF_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))

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

    if (on_root) then

      allocate (Lo_qb1_q_qb2(num_bands, num_bands))
      allocate (L_qb1_q_qb2(num_wann, num_wann))
      allocate (FF_q(num_wann, num_wann, num_kpts, 3, 3))

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

      open (newunit=uIu_in, file=TRIM(seedname)//".uIu", form='unformatted', &
            status='old', action='read', err=107)
      write (stdout, '(/a)', advance='no') &
        ' Reading uIu overlaps from '//trim(seedname)//'.uIu in get_FF_R: '
      read (uIu_in, err=108, end=108) header
      write (stdout, '(a)') trim(header)
      read (uIu_in, err=108, end=108) nb_tmp, nkp_tmp, nntot_tmp
      if (nb_tmp .ne. num_bands) then
        call set_error_fatal(error, trim(seedname)//'.uIu has not the right number of bands', comm)
        return
      end if
      if (nkp_tmp .ne. num_kpts) then
        call set_error_fatal(error, trim(seedname)//'.uIu has not the right number of k-points', comm)
        return
      end if
      if (nntot_tmp .ne. kmesh_info%nntot) then
        call set_error_fatal(error, trim(seedname)//'.uIu has not the right number of nearest neighbours', comm)
        return
      end if

      FF_q = cmplx_0
      do ik = 1, num_kpts
        do nn2 = 1, kmesh_info%nntot
          qb2 = kmesh_info%nnlist(ik, nn2)

          call get_win_min(num_bands, dis_manifold, qb2, winmin_qb2, have_disentangled)
          do nn1 = 1, kmesh_info%nntot
            qb1 = kmesh_info%nnlist(ik, nn1)
            call get_win_min(num_bands, dis_manifold, qb1, winmin_qb1, have_disentangled)
            !
            ! Read from .uIu file the matrices <u_{q+b1}|u_{q+b2}>
            ! between the original ab initio eigenstates
            !
            !               do m=1,num_bands
            !                  do n=1,num_bands
            !                     read(uIu_in,err=108,end=108) Lo_qb1_q_qb2(m,n)
            !                  end do
            !               end do
            read (uIu_in, err=108, end=108) &
              ((Lo_qb1_q_qb2(n, m), n=1, num_bands), m=1, num_bands)
            !
            ! **************************************************************
            ! 2013-08-09: Do we need to take a transpose here?! SEE get_CC_R
            Lo_qb1_q_qb2 = transpose(Lo_qb1_q_qb2) ! added 2013-08-09 (?)
            ! **************************************************************
            !
            ! Transform to projected subspace, Wannier gauge
            !
            L_qb1_q_qb2(:, :) = cmplx_0
            do m = 1, num_wann
              do n = 1, num_wann
                do i = 1, num_states(qb1)
                  ii = winmin_qb1 + i - 1
                  do j = 1, num_states(qb2)
                    jj = winmin_qb2 + j - 1
                    L_qb1_q_qb2(n, m) = L_qb1_q_qb2(n, m) &
                                        + conjg(v_matrix(i, n, qb1)) &
                                        *Lo_qb1_q_qb2(ii, jj) &
                                        *v_matrix(j, m, qb2)
                  end do
                end do
              end do
            end do
            do b = 1, 3
              do a = 1, b
                FF_q(:, :, ik, a, b) = FF_q(:, :, ik, a, b) &
                                       + kmesh_info%wb(nn1)*kmesh_info%bk(a, nn1, ik) &
                                       *kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*L_qb1_q_qb2(:, :)
              end do
            end do
          end do !nn1
        end do !nn2
        do b = 1, 3
          do a = 1, b
            FF_q(:, :, ik, b, a) = conjg(transpose(FF_q(:, :, ik, a, b)))
          end do
        end do
      end do !ik

      close (uIu_in)

      do b = 1, 3
        do a = 1, 3
          call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
                              FF_q(:, :, :, a, b), FF_R_temp(:, :, :, a, b))
        end do
      end do

      do b = 1, 3
        do a = 1, 3
          call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
                                     FF_R_temp(:, :, :, a, b), FF_R(:, :, :, a, b))
        end do
      end do

    end if !on_root

    call comms_bcast(FF_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
    if (allocated(error)) return

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

    deallocate (FF_R_temp)

107 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.uIu', comm)
    return
108 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.uIu', comm)
    return

  end subroutine get_FF_R