get_AA_R Subroutine

public subroutine get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wann_data, AA_R, v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, timer, error, comm)

Uses

  • proc~~get_aa_r~~UsesGraph proc~get_aa_r get_AA_R module~w90_postw90_types w90_postw90_types proc~get_aa_r->module~w90_postw90_types module~w90_types w90_types proc~get_aa_r->module~w90_types module~w90_comms w90_comms module~w90_postw90_types->module~w90_comms module~w90_constants w90_constants module~w90_postw90_types->module~w90_constants module~w90_types->module~w90_constants module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base

AA_a(R) = <0|r_a|R> is the Fourier transform of the Berrry connection AA_a(k) = i (a=x,y,z)

Arguments

Type IntentOptional Attributes Name
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
type(kmesh_info_type), intent(in) :: kmesh_info
real(kind=dp), intent(in) :: kpt_latt(:,:)
type(print_output_type), intent(in) :: print_output
type(wannier_data_type), intent(in) :: wann_data
complex(kind=dp), intent(inout), allocatable :: AA_R(:,:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
real(kind=dp), intent(in) :: eigval(:,:)
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(ws_distance_type), intent(in) :: ws_distance
type(ws_region_type), intent(in) :: ws_region
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~get_aa_r~~CallsGraph proc~get_aa_r get_AA_R interface~comms_bcast comms_bcast proc~get_aa_r->interface~comms_bcast interface~comms_reduce comms_reduce proc~get_aa_r->interface~comms_reduce interface~comms_scatterv comms_scatterv proc~get_aa_r->interface~comms_scatterv proc~comms_array_split comms_array_split proc~get_aa_r->proc~comms_array_split proc~fourier_loc_q_to_r fourier_loc_q_to_R proc~get_aa_r->proc~fourier_loc_q_to_r proc~get_gauge_overlap_matrix get_gauge_overlap_matrix proc~get_aa_r->proc~get_gauge_overlap_matrix proc~io_stopwatch_start io_stopwatch_start proc~get_aa_r->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~get_aa_r->proc~io_stopwatch_stop proc~mpirank mpirank proc~get_aa_r->proc~mpirank proc~mpisize mpisize proc~get_aa_r->proc~mpisize proc~operator_wigner_setup operator_wigner_setup proc~get_aa_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_aa_r->proc~set_error_fatal proc~set_error_file set_error_file proc~get_aa_r->proc~set_error_file proc~comms_bcast_char comms_bcast_char interface~comms_bcast->proc~comms_bcast_char proc~comms_bcast_cmplx comms_bcast_cmplx interface~comms_bcast->proc~comms_bcast_cmplx proc~comms_bcast_int comms_bcast_int interface~comms_bcast->proc~comms_bcast_int proc~comms_bcast_logical comms_bcast_logical interface~comms_bcast->proc~comms_bcast_logical proc~comms_bcast_real comms_bcast_real interface~comms_bcast->proc~comms_bcast_real proc~comms_reduce_cmplx comms_reduce_cmplx interface~comms_reduce->proc~comms_reduce_cmplx proc~comms_reduce_int comms_reduce_int interface~comms_reduce->proc~comms_reduce_int proc~comms_reduce_real comms_reduce_real interface~comms_reduce->proc~comms_reduce_real proc~comms_scatterv_cmplx_3 comms_scatterv_cmplx_3 interface~comms_scatterv->proc~comms_scatterv_cmplx_3 proc~comms_scatterv_cmplx_4 comms_scatterv_cmplx_4 interface~comms_scatterv->proc~comms_scatterv_cmplx_4 proc~comms_scatterv_int_1 comms_scatterv_int_1 interface~comms_scatterv->proc~comms_scatterv_int_1 proc~comms_scatterv_int_2 comms_scatterv_int_2 interface~comms_scatterv->proc~comms_scatterv_int_2 proc~comms_scatterv_int_3 comms_scatterv_int_3 interface~comms_scatterv->proc~comms_scatterv_int_3 proc~comms_scatterv_real_1 comms_scatterv_real_1 interface~comms_scatterv->proc~comms_scatterv_real_1 proc~comms_scatterv_real_2 comms_scatterv_real_2 interface~comms_scatterv->proc~comms_scatterv_real_2 proc~comms_scatterv_real_3 comms_scatterv_real_3 interface~comms_scatterv->proc~comms_scatterv_real_3 proc~comms_array_split->proc~mpisize proc~get_win_min get_win_min proc~get_gauge_overlap_matrix->proc~get_win_min proc~utility_zgemmm utility_zgemmm proc~get_gauge_overlap_matrix->proc~utility_zgemmm proc~comms_sync_error comms_sync_error proc~set_error_fatal->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_fatal->proc~set_base_error proc~set_error_file->proc~comms_sync_error proc~set_error_file->proc~set_base_error proc~comms_bcast_char->proc~comms_sync_error proc~comms_no_sync_bcast_char comms_no_sync_bcast_char proc~comms_bcast_char->proc~comms_no_sync_bcast_char proc~comms_bcast_cmplx->proc~comms_sync_error proc~comms_no_sync_bcast_cmplx comms_no_sync_bcast_cmplx proc~comms_bcast_cmplx->proc~comms_no_sync_bcast_cmplx proc~comms_bcast_int->proc~comms_sync_error proc~comms_no_sync_bcast_int comms_no_sync_bcast_int proc~comms_bcast_int->proc~comms_no_sync_bcast_int proc~comms_bcast_logical->proc~comms_sync_error proc~comms_no_sync_bcast_logical comms_no_sync_bcast_logical proc~comms_bcast_logical->proc~comms_no_sync_bcast_logical proc~comms_bcast_real->proc~comms_sync_error proc~comms_no_sync_bcast_real comms_no_sync_bcast_real proc~comms_bcast_real->proc~comms_no_sync_bcast_real proc~comms_reduce_cmplx->proc~comms_sync_error proc~comms_no_sync_reduce_cmplx comms_no_sync_reduce_cmplx proc~comms_reduce_cmplx->proc~comms_no_sync_reduce_cmplx proc~comms_reduce_int->proc~comms_sync_error proc~comms_no_sync_reduce_int comms_no_sync_reduce_int proc~comms_reduce_int->proc~comms_no_sync_reduce_int proc~comms_reduce_real->proc~comms_sync_error proc~comms_no_sync_reduce_real comms_no_sync_reduce_real proc~comms_reduce_real->proc~comms_no_sync_reduce_real proc~comms_scatterv_cmplx_3->proc~comms_sync_error proc~comms_no_sync_scatterv_cmplx_3 comms_no_sync_scatterv_cmplx_3 proc~comms_scatterv_cmplx_3->proc~comms_no_sync_scatterv_cmplx_3 proc~comms_scatterv_cmplx_4->proc~comms_sync_error proc~comms_no_sync_scatterv_cmplx_4 comms_no_sync_scatterv_cmplx_4 proc~comms_scatterv_cmplx_4->proc~comms_no_sync_scatterv_cmplx_4 proc~comms_scatterv_int_1->proc~comms_sync_error proc~comms_no_sync_scatterv_int_1 comms_no_sync_scatterv_int_1 proc~comms_scatterv_int_1->proc~comms_no_sync_scatterv_int_1 proc~comms_scatterv_int_2->proc~comms_sync_error proc~comms_no_sync_scatterv_int_2 comms_no_sync_scatterv_int_2 proc~comms_scatterv_int_2->proc~comms_no_sync_scatterv_int_2 proc~comms_scatterv_int_3->proc~comms_sync_error proc~comms_no_sync_scatterv_int_3 comms_no_sync_scatterv_int_3 proc~comms_scatterv_int_3->proc~comms_no_sync_scatterv_int_3 proc~comms_scatterv_real_1->proc~comms_sync_error proc~comms_no_sync_scatterv_real_1 comms_no_sync_scatterv_real_1 proc~comms_scatterv_real_1->proc~comms_no_sync_scatterv_real_1 proc~comms_scatterv_real_2->proc~comms_sync_error proc~comms_no_sync_scatterv_real_2 comms_no_sync_scatterv_real_2 proc~comms_scatterv_real_2->proc~comms_no_sync_scatterv_real_2 proc~comms_scatterv_real_3->proc~comms_sync_error proc~comms_no_sync_scatterv_real_3 comms_no_sync_scatterv_real_3 proc~comms_scatterv_real_3->proc~comms_no_sync_scatterv_real_3 proc~utility_zgemm_new utility_zgemm_new proc~utility_zgemmm->proc~utility_zgemm_new zgemm zgemm proc~utility_zgemm_new->zgemm

Called by

proc~~get_aa_r~~CalledByGraph proc~get_aa_r get_AA_R proc~berry_main berry_main proc~berry_main->proc~get_aa_r proc~berry_get_sc_klist berry_get_sc_klist proc~berry_main->proc~berry_get_sc_klist proc~gyrotropic_main gyrotropic_main proc~gyrotropic_main->proc~get_aa_r proc~k_path k_path proc~k_path->proc~get_aa_r proc~k_slice k_slice proc~k_slice->proc~get_aa_r 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_aa_r proc~berry_get_sc_klist->proc~wham_get_eig_uu_hh_aa_sc_tb_conv program~postw90 postw90 program~postw90->proc~berry_main program~postw90->proc~gyrotropic_main program~postw90->proc~k_path program~postw90->proc~k_slice

Source Code

  subroutine get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wann_data, AA_R, &
                      v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
                      num_wann, have_disentangled, seedname, stdout, timer, error, comm)
    !================================================
    !
    !! AA_a(R) = <0|r_a|R> is the Fourier transform
    !! of the Berrry connection AA_a(k) = i<u|del_a u>
    !! (a=x,y,z)
    !
    !================================================

    use w90_postw90_types, only: pw90_berry_mod_type, pw90_oper_read_type, pw90_spin_hall_type, &
                                 wigner_seitz_type
    use w90_types, only: dis_manifold_type, kmesh_info_type, print_output_type, timer_list_type, &
                         ws_distance_type, ws_region_type, wannier_data_type

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(in) :: pw90_berry
    type(dis_manifold_type), intent(in)   :: dis_manifold
    type(kmesh_info_type), intent(in)     :: kmesh_info
    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(print_output_type), intent(in)   :: print_output
    type(wannier_data_type), intent(in) :: wann_data
    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) :: eigval(:, :)
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    complex(kind=dp), intent(in) :: v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
    complex(kind=dp), allocatable :: AA_R_temp(:, :, :)
    complex(kind=dp), allocatable :: AA_R_b(:, :, :, :)

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

    ! local variables
    complex(kind=dp), allocatable :: AA_q_b(:, :, :, :, :)
    complex(kind=dp), allocatable :: AA_q(:, :, :, :)
    complex(kind=dp), allocatable :: AA_q_b_diag(:, :, :)
    complex(kind=dp), allocatable :: AA_q_loc(:, :, :)
    complex(kind=dp), allocatable :: S_o(:, :)
    complex(kind=dp), allocatable :: S(:, :)
    integer                       :: n, m, i, j, &
                                     ik, ik2, ik_prev, nn, inn, nnl, nnm, nnn, &
                                     idir, ncount, nn_count, mmn_in, &
                                     nb_tmp, nkp_tmp, nntot_tmp, file_unit, &
                                     ir, io, w
    integer, allocatable          :: num_states(:)
    integer, allocatable          :: counts(:), displs(:)
    real(kind=dp)                 :: m_real, m_imag, rdum1_real, rdum1_imag, &
                                     rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
    real(kind=dp), allocatable    :: r0(:, :, :)
    complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
    logical                       :: nn_found
    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_AA_R', timer)

    if (.not. allocated(wigner_seitz%wannier_centres_from_AA_R)) then
      allocate (wigner_seitz%wannier_centres_from_AA_R(3, num_wann))
    end if

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

    ! Real-space position matrix elements calculated by Fourier
    ! transforming overlap matrices defined on the ab-initio
    ! reciprocal mesh
    !
    ! Do everything on root, broadcast AA_R at the end (smaller than S_o)
    !
    if (on_root) then
      allocate (AA_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, 3))
      AA_R = cmplx_0
    else
      allocate (AA_q_b(1, 1, 1, kmesh_info%nntot, 3))
    end if
    !
    if (on_root) then
      allocate (S_o(num_bands, num_bands))
      allocate (S(num_wann, num_wann))
      allocate (AA_q_b_diag(num_wann, kmesh_info%nntot, 3))

      allocate (num_states(num_kpts))

      wigner_seitz%wannier_centres_from_AA_R(:, :) = 0.d0

      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=mmn_in, file=trim(seedname)//'.mmn', &
            form='formatted', status='old', action='read', err=101)
      write (stdout, '(/a)', advance='no') &
        ' Reading overlaps from '//trim(seedname)//'.mmn in get_AA_R   : '
      ! Read the comment line (header)
      read (mmn_in, '(a)', err=102, end=102) header
      write (stdout, '(a)') trim(header)
      ! Read the number of bands, k-points and nearest neighbours
      read (mmn_in, *, err=102, end=102) nb_tmp, nkp_tmp, nntot_tmp
      ! Checks
      if (nb_tmp .ne. num_bands) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of bands', comm)
        return
      end if
      if (nkp_tmp .ne. num_kpts) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of k-points', comm)
        return
      end if
      if (nntot_tmp .ne. kmesh_info%nntot) then
        call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of nearest neighbours', comm)
        return
      end if

      AA_q_b = cmplx_0
      ik_prev = 0

      ! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
      do ncount = 1, num_kpts*kmesh_info%nntot
        !
        !Read from .mmn file the original overlap matrix
        ! S_o=<u_ik|u_ik2> between ab initio eigenstates
        !
        read (mmn_in, *, err=102, end=102) ik, ik2, nnl, nnm, nnn
        do n = 1, num_bands
          do m = 1, num_bands
            read (mmn_in, *, err=102, end=102) m_real, m_imag
            S_o(m, n) = cmplx(m_real, m_imag, kind=dp)
          end do
        end do

        !debug
        !OK
        !if(ik.ne.ik_prev .and.ik_prev.ne.0) then
        !   if(nn_count.ne.nntot)&
        !        write(stdout,*) 'something wrong in get_AA_R!'
        !endif
        !enddebug
        if (ik .ne. ik_prev) nn_count = 0
        nn = 0
        nn_found = .false.
        do inn = 1, kmesh_info%nntot
          if ((ik2 .eq. kmesh_info%nnlist(ik, inn)) .and. &
              (nnl .eq. kmesh_info%nncell(1, ik, inn)) .and. &
              (nnm .eq. kmesh_info%nncell(2, ik, inn)) .and. &
              (nnn .eq. kmesh_info%nncell(3, ik, inn))) then
            if (.not. nn_found) then
              nn_found = .true.
              nn = inn
            else
              call set_error_fatal(error, 'Error reading '//trim(seedname)//'.mmn.&
                   & More than one matching nearest neighbour found', comm)
              return
            end if
          end if
        end do
        if (nn .eq. 0) then
          write (stdout, '(/a,i8,2i5,i4,2x,3i3)') ' Error reading '//trim(seedname)//'.mmn:', &
            ncount, ik, ik2, nn, nnl, nnm, nnn
          call set_error_fatal(error, 'Neighbour not found', comm)
          return
        end if
        nn_count = nn_count + 1 !Check: can also be place after nn=inn (?)

        ! Wannier-gauge overlap matrix S in the projected subspace
        !
        call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                      ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
                                      num_states(kmesh_info%nnlist(ik, nn)), S_o, &
                                      have_disentangled, S)

        ! save the wannier centers (diagonals of AA_R_temp) to wannier_centres_from_AA_R
        ! used in pw90common_fourier_R_to_k_new_second_d_TB_conv
        do i = 1, num_wann
          wigner_seitz%wannier_centres_from_AA_R(:, i) = &
            wigner_seitz%wannier_centres_from_AA_R(:, i) &
            - kmesh_info%wb(nn)*kmesh_info%bk(:, nn, ik)*aimag(log(S(i, i)))/num_kpts
        end do

        ! Berry connection matrix
        !
        if (pw90_berry%transl_inv .and. ik .ne. ik_prev) AA_q_b_diag(:, :, :) = cmplx_0

        nno = nn
        if (pw90_berry%transl_inv_full) nno = kmesh_info%nninv(nn, ik) ! reorder AA_q_b nn indices required for transl_inv_full method

        do idir = 1, 3
          AA_q_b(:, :, ik, nno, idir) = AA_q_b(:, :, ik, nno, idir) &
                                        + cmplx_i*kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik)*S(:, :)
          if (pw90_berry%transl_inv) then
            !
            ! Rewrite band-diagonal elements a la Eq.(31) of MV97
            !
            do i = 1, num_wann
              AA_q_b_diag(i, nno, idir) = AA_q_b_diag(i, nno, idir) &
                                          - kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
                                          *aimag(log(S(i, i)))
            end do
          end if
        end do

        do idir = 1, 3
          if (pw90_berry%transl_inv) then
            do n = 1, num_wann
              AA_q_b(n, n, ik, nno, idir) = AA_q_b_diag(n, nno, idir)
            end do
          end if
        end do

        ik_prev = ik
      end do !ncount

      close (mmn_in)

      if (sum((wigner_seitz%wannier_centres_from_AA_R - wann_data%centres)**2) > 1.0e-08) then
        if (pw90_berry%guiding_centres) then
          write (stdout, '(/a)', advance='no') &
            ' Computed and read Wannier centres different. This can happen for guiding_centres=T'
          wigner_seitz%wannier_centres_from_AA_R = wann_data%centres
        else
          call set_error_fatal(error, 'Computed and read Wannier centres different.', comm)
        end if
      end if

      if (pw90_berry%transl_inv_full) then
        allocate (r0(num_wann, num_wann, 3))
        allocate (phase1(num_wann, num_wann))
        do j = 1, num_wann
          do i = 1, num_wann
            r0(i, j, :) = (wigner_seitz%wannier_centres_from_AA_R(:, i) + &
                           wigner_seitz%wannier_centres_from_AA_R(:, j))/2.0_dp
          end do
        end do

        do nn = 1, kmesh_info%nntot
          do ik = 1, num_kpts
            phase1 = (r0(:, :, 1)*kmesh_info%bk(1, nn, ik) + &
                      r0(:, :, 2)*kmesh_info%bk(2, nn, ik) + &
                      r0(:, :, 3)*kmesh_info%bk(3, nn, ik))
            phase1 = exp(cmplx_i*phase1)

            nno = kmesh_info%nninv(nn, ik)
            do idir = 1, 3
              AA_q_b(:, :, ik, nno, idir) = AA_q_b(:, :, ik, nno, idir)*phase1(:, :)
            end do
          end do ! ik
        end do ! nn
        deallocate (phase1)
      end if

    end if !on_root

    if (pw90_berry%transl_inv_full) then
      allocate (counts(0:mpisize(comm) - 1))
      allocate (displs(0:mpisize(comm) - 1))

      w = num_wann*num_wann
      call comms_array_split(num_kpts, counts, displs, comm)
      allocate (AA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
      allocate (AA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
      if (on_root) then
        allocate (AA_R_b(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
        allocate (phase2(wigner_seitz%nrpts_pw90))
      end if

      do idir = 1, 3
        do nn = 1, kmesh_info%nntot
          call comms_scatterv(AA_q_loc, w*counts(mpirank(comm)), AA_q_b(:, :, :, nn, idir), w*counts, w*displs, error, comm)
          call fourier_loc_q_to_R(num_kpts, counts, displs, mpirank(comm), wigner_seitz%nrpts, wigner_seitz%irvec, &
                                  kpt_latt, AA_q_loc, AA_R_temp)
          call comms_reduce(AA_R_temp(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, 'SUM', error, comm)

          if (on_root) then
            ! Apply degeneracy factor and reorder according to the wigner-seitz vectors
            call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, AA_R_temp, AA_R_b(:, :, :, idir))

            phase2 = -0.5_dp*(wigner_seitz%crvec_pw90(1, :)*kmesh_info%bk(1, nn, 1) + &
                              wigner_seitz%crvec_pw90(2, :)*kmesh_info%bk(2, nn, 1) + &
                              wigner_seitz%crvec_pw90(3, :)*kmesh_info%bk(3, nn, 1))

            phase2 = exp(cmplx_i*phase2)
            AA_R(:, :, :, idir) = AA_R(:, :, :, idir) + AA_R_b(:, :, :, idir)*spread(spread(phase2, 1, num_wann), 1, num_wann)
          end if
        end do
      end do

      deallocate (AA_q_loc)
      deallocate (AA_R_temp)

      if (on_root) then
        deallocate (phase2)
        deallocate (AA_q_b)
        deallocate (AA_R_b)

        do ir = 1, wigner_seitz%nrpts_pw90
          if ((wigner_seitz%irvec_pw90(1, ir) .eq. 0) .and. &
              (wigner_seitz%irvec_pw90(2, ir) .eq. 0) .and. &
              (wigner_seitz%irvec_pw90(3, ir) .eq. 0)) then
            do i = 1, num_wann
              AA_R(i, i, ir, :) = wigner_seitz%wannier_centres_from_AA_R(:, i)
            end do
            exit
          end if
        end do
      end if
    else
      allocate (AA_q(num_wann, num_wann, num_kpts, 3))

      if (on_root) then
        AA_q = sum(AA_q_b, 4)
        deallocate (AA_q_b)

        ! Since Eq.(44) WYSV06 does not preserve the Hermiticity of the
        ! Berry potential matrix, take Hermitean part (whether this
        ! makes a difference or not for e.g. the AHC, depends on which
        ! expression is used to evaluate the Berry curvature.
        ! See comments in berry_wanint.F90)
        !
        do idir = 1, 3
          do ik = 1, num_kpts
            AA_q(:, :, ik, idir) = &
              0.5_dp*(AA_q(:, :, ik, idir) &
                      + conjg(transpose(AA_q(:, :, ik, idir))))
          end do
        end do
      end if
      !
      allocate (counts(0:mpisize(comm) - 1))
      allocate (displs(0:mpisize(comm) - 1))

      w = num_wann*num_wann
      call comms_array_split(num_kpts, counts, displs, comm)
      allocate (AA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
      allocate (AA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))

      do idir = 1, 3
        call comms_scatterv(AA_q_loc, w*counts(mpirank(comm)), AA_q(:, :, :, idir), w*counts, w*displs, error, comm)
        call fourier_loc_q_to_R(num_kpts, counts, displs, mpirank(comm), wigner_seitz%nrpts, wigner_seitz%irvec, &
                                kpt_latt, AA_q_loc, AA_R_temp)
        call comms_reduce(AA_R_temp(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, 'SUM', error, comm)

        if (on_root) then
          ! Apply degeneracy factor and reorder according to the wigner-seitz vectors
          call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, AA_R_temp, AA_R(:, :, :, idir))
        end if
      end do

      deallocate (AA_q_loc)
      deallocate (AA_q)
      deallocate (AA_R_temp)
    end if

    call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3, error, comm)
    call comms_bcast(wigner_seitz%wannier_centres_from_AA_R(1, 1), num_wann*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_AA_R', timer)
    return

101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
    return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
    return
103 call set_error_file(error, 'Error in get_AA_R: problem opening file '//trim(seedname)//'_AA_R.dat', comm)
    return

  end subroutine get_AA_R