get_CC_R Subroutine

public subroutine get_CC_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, HH_R, BB_R, CC_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, timer, error, comm)

Uses

  • proc~~get_cc_r~~UsesGraph proc~get_cc_r get_CC_R module~w90_postw90_types w90_postw90_types proc~get_cc_r->module~w90_postw90_types module~w90_types w90_types proc~get_cc_r->module~w90_types module~w90_utility w90_utility proc~get_cc_r->module~w90_utility 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_utility->module~w90_comms module~w90_utility->module~w90_constants module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base

CC_ab(R) = <0|r_a.H.(r-R)_b|R> is the Fourier transform of CC_ab(k) = (a,b=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(pw90_oper_read_type), intent(in) :: pw90_oper_read
complex(kind=dp), intent(inout), allocatable :: HH_R(:,:,:)
complex(kind=dp), intent(inout), allocatable :: BB_R(:,:,:,:)
complex(kind=dp), intent(inout), allocatable :: CC_R(:,:,:,:,:)
complex(kind=dp), intent(in) :: v_matrix(:,:,:)
real(kind=dp), intent(in) :: eigval(:,:)
real(kind=dp), intent(in) :: scissors_shift
type(wigner_seitz_type), intent(in) :: 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_cc_r~~CallsGraph proc~get_cc_r get_CC_R interface~comms_bcast comms_bcast proc~get_cc_r->interface~comms_bcast interface~comms_reduce comms_reduce proc~get_cc_r->interface~comms_reduce interface~comms_scatterv comms_scatterv proc~get_cc_r->interface~comms_scatterv proc~comms_array_split comms_array_split proc~get_cc_r->proc~comms_array_split proc~fourier_loc_q_to_r fourier_loc_q_to_R proc~get_cc_r->proc~fourier_loc_q_to_r proc~get_gauge_overlap_matrix get_gauge_overlap_matrix proc~get_cc_r->proc~get_gauge_overlap_matrix proc~get_win_min get_win_min proc~get_cc_r->proc~get_win_min proc~io_stopwatch_start io_stopwatch_start proc~get_cc_r->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~get_cc_r->proc~io_stopwatch_stop proc~mpirank mpirank proc~get_cc_r->proc~mpirank proc~mpisize mpisize proc~get_cc_r->proc~mpisize proc~operator_wigner_setup operator_wigner_setup proc~get_cc_r->proc~operator_wigner_setup proc~set_error_fatal set_error_fatal proc~get_cc_r->proc~set_error_fatal proc~set_error_file set_error_file proc~get_cc_r->proc~set_error_file proc~utility_compar utility_compar proc~get_cc_r->proc~utility_compar 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_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_cc_r~~CalledByGraph proc~get_cc_r get_CC_R proc~berry_main berry_main proc~berry_main->proc~get_cc_r proc~gyrotropic_main gyrotropic_main proc~gyrotropic_main->proc~get_cc_r proc~k_path k_path proc~k_path->proc~get_cc_r proc~k_slice k_slice proc~k_slice->proc~get_cc_r 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_CC_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, HH_R, &
                      BB_R, CC_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, &
                      num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, timer, error, comm)
    !================================================
    !
    !! CC_ab(R) = <0|r_a.H.(r-R)_b|R> is the Fourier transform of
    !! CC_ab(k) = <del_a u|H|del_b u> (a,b=x,y,z)
    !
    !================================================

    use w90_postw90_types, only: pw90_berry_mod_type, pw90_oper_read_type, 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
    use w90_utility, only: utility_compar

    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(pw90_oper_read_type), intent(in) :: pw90_oper_read
    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) :: eigval(:, :)
    real(kind=dp), intent(in) :: scissors_shift
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    complex(kind=dp), intent(in) :: v_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :)
    complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :) ! <0|r_alpha.H(r-R)_beta|R>
    complex(kind=dp), allocatable :: CC_R_temp(:, :, :)
    complex(kind=dp), allocatable :: CC_R_b(:, :, :, :, :)

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

    ! local variables
    integer          :: m, n, a, b, nn1, nn2, ik, nb_tmp, nkp_tmp, i, j, ir, ir2, &
                        nntot_tmp, uHu_in, qb1, qb2, winmin_qb1, winmin_qb2, &
                        ifpos, ifneg, w

    integer, allocatable          :: num_states(:)
    integer, allocatable          :: counts(:), displs(:)
    complex(kind=dp), allocatable :: CC_q(:, :, :, :, :)
    complex(kind=dp), allocatable :: CC_q_b(:, :, :, :, :, :, :)
    complex(kind=dp), allocatable :: CC_q_loc(:, :, :)
    complex(kind=dp), allocatable :: Ho_qb1_q_qb2(:, :)
    complex(kind=dp), allocatable :: H_qb1_q_qb2(:, :)
    real(kind=dp), allocatable    :: r0(:, :, :)
    complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
    real(kind=dp)                 :: c_real, c_img
    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_CC_R', timer)

    if (.not. allocated(CC_R)) then
      allocate (CC_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_CC_R', timer)
      return
    end if

    if (on_root) then
      allocate (CC_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, kmesh_info%nntot, 3, 3))
      CC_R = cmplx_0
    else
      allocate (CC_q_b(1, 1, 1, kmesh_info%nntot, kmesh_info%nntot, 3, 3))
    end if

    if (on_root) then

      if (abs(scissors_shift) > 1.0e-7_dp) then
        call set_error_fatal(error, 'Error: scissors correction not yet implemented for CC_R', comm)
        return
      end if

      allocate (Ho_qb1_q_qb2(num_bands, num_bands))
      allocate (H_qb1_q_qb2(num_wann, num_wann))

      allocate (phase1(num_wann, num_wann))
      if (pw90_berry%transl_inv_full) then
        allocate (r0(num_wann, num_wann, 3))
        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
      end if

      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

      if (pw90_oper_read%uHu_formatted) then
        open (newunit=uHu_in, file=trim(seedname)//".uHu", form='formatted', &
              status='old', action='read', err=105)
        write (stdout, '(/a)', advance='no') &
          ' Reading uHu overlaps from '//trim(seedname)//'.uHu in get_CC_R: '
        read (uHu_in, *, err=106, end=106) header
        write (stdout, '(a)') trim(header)
        read (uHu_in, *, err=106, end=106) nb_tmp, nkp_tmp, nntot_tmp
      else
        open (newunit=uHu_in, file=trim(seedname)//".uHu", form='unformatted', &
              status='old', action='read', err=105)
        write (stdout, '(/a)', advance='no') &
          ' Reading uHu overlaps from '//trim(seedname)//'.uHu in get_CC_R: '
        read (uHu_in, err=106, end=106) header
        write (stdout, '(a)') trim(header)
        read (uHu_in, err=106, end=106) nb_tmp, nkp_tmp, nntot_tmp
      end if
      if (nb_tmp .ne. num_bands) then
        call set_error_fatal(error, trim(seedname)//'.uHu 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)//'.uHu 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)//'.uHu has not the right number of nearest neighbours', comm)
        return
      end if

      CC_q_b = 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 .uHu file the matrices <u_{q+b1}|H_q|u_{q+b2}>
            ! between the original ab initio eigenstates
            !
            if (pw90_oper_read%uHu_formatted) then
              do m = 1, num_bands
                do n = 1, num_bands
                  read (uHu_in, *, err=106, end=106) c_real, c_img
                  Ho_qb1_q_qb2(n, m) = cmplx(c_real, c_img, dp)
                end do
              end do
            else
              read (uHu_in, err=106, end=106) &
                ((Ho_qb1_q_qb2(n, m), n=1, num_bands), m=1, num_bands)
            end if
            ! pw2wannier90 is coded a bit strangely, so here we take the transpose
            Ho_qb1_q_qb2 = transpose(Ho_qb1_q_qb2)
            ! old code here
            !do m=1,num_bands
            !   do n=1,num_bands
            !      read(uHu_in,err=106,end=106) Ho_qb1_q_qb2(m,n)
            !   end do
            !end do
            !
            ! Transform to projected subspace, Wannier gauge
            !

            call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
                                          qb1, num_states(qb1), qb2, num_states(qb2), &
                                          Ho_qb1_q_qb2, have_disentangled, H_qb1_q_qb2)

            if (pw90_berry%transl_inv_full) then
              phase1 = -(r0(:, :, 1)*kmesh_info%bk(1, nn1, ik) + &
                         r0(:, :, 2)*kmesh_info%bk(2, nn1, ik) + &
                         r0(:, :, 3)*kmesh_info%bk(3, nn1, ik)) &
                       + (r0(:, :, 1)*kmesh_info%bk(1, nn2, ik) + &
                          r0(:, :, 2)*kmesh_info%bk(2, nn2, ik) + &
                          r0(:, :, 3)*kmesh_info%bk(3, nn2, ik))

              phase1 = exp(cmplx_i*phase1)
            else
              phase1 = cmplx_1
            end if

            do b = 1, 3
              do a = 1, b
                nn1o = kmesh_info%nninv(nn1, ik)
                nn2o = kmesh_info%nninv(nn2, ik)
                CC_q_b(:, :, ik, nn1o, nn2o, a, b) = CC_q_b(:, :, ik, nn1o, nn2o, a, b) &
                                                     + phase1(:, :)*kmesh_info%wb(nn1)*kmesh_info%bk(a, nn1, ik) &
                                                     *kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*H_qb1_q_qb2(:, :)
              end do
            end do

          end do !nn1
        end do !nn2
      end do !ik

      close (uHu_in)

      deallocate (phase1)

    end if !on_root

    if (pw90_berry%transl_inv_full) then
      if (.not. allocated(HH_R)) then
        call set_error_fatal(error, 'transl_inv_full=T for CC_R needs HH_R', comm)
      end if

      if (.not. allocated(BB_R)) then
        call set_error_fatal(error, 'transl_inv_full=T for CC_R needs BB_R', comm)
      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 (CC_q_loc(num_wann, num_wann, counts(mpirank(comm))))
      allocate (CC_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
      if (on_root) then
        allocate (CC_R_b(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3, 3))
        allocate (phase2(wigner_seitz%nrpts_pw90))
      end if

      do b = 1, 3
        do a = 1, 3
          do nn2 = 1, kmesh_info%nntot
            do nn1 = 1, kmesh_info%nntot
              call comms_scatterv(CC_q_loc, w*counts(mpirank(comm)), CC_q_b(:, :, :, nn1, nn2, a, b), &
                                  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, CC_q_loc, CC_R_temp)
              call comms_reduce(CC_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, &
                                           CC_R_temp, CC_R_b(:, :, :, a, b))

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

                phase2 = exp(cmplx_i*phase2)
                CC_R(:, :, :, a, b) = CC_R(:, :, :, a, b) + CC_R_b(:, :, :, a, b)* &
                                      spread(spread(phase2, 1, num_wann), 1, num_wann)
              end if
            end do
          end do
        end do
      end do

      deallocate (CC_q_loc)
      deallocate (CC_R_temp)

      if (on_root) then
        deallocate (phase2)
        deallocate (CC_q_b)
        deallocate (CC_R_b)

        do b = 1, 3
          do a = 1, 3
            do ir = 1, wigner_seitz%nrpts_pw90
              CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
                                     (r0(:, :, a) + 0.5_dp*wigner_seitz%crvec_pw90(a, ir))* &
                                     BB_R(:, :, ir, b)
              do ir2 = 1, wigner_seitz%nrpts_pw90
                call utility_compar(wigner_seitz%crvec_pw90(1, ir), &
                                    wigner_seitz%crvec_pw90(1, ir2), ifpos, ifneg)
                if (ifneg .eq. 1) then
                  CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
                                         conjg(transpose(BB_R(:, :, ir2, a)))* &
                                         (r0(:, :, b) - 0.5_dp*wigner_seitz%crvec_pw90(b, ir))
                  exit
                end if
              end do
              CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
                                     (r0(:, :, a) + 0.5_dp*wigner_seitz%crvec_pw90(a, ir))* &
                                     wigner_seitz%crvec_pw90(b, ir)*HH_R(:, :, ir)
            end do
          end do
        end do
      end if
    else
      allocate (CC_q(num_wann, num_wann, num_kpts, 3, 3))

      if (on_root) then
        CC_q = sum(sum(CC_q_b, 5), 4)
        deallocate (CC_q_b)
        !
        do b = 1, 3
          do a = 1, b
            do ik = 1, num_kpts
              CC_q(:, :, ik, b, a) = conjg(transpose(CC_q(:, :, ik, a, b)))
            end do
          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 (CC_q_loc(num_wann, num_wann, counts(mpirank(comm))))
      allocate (CC_R_temp(num_wann, num_wann, wigner_seitz%nrpts))

      do b = 1, 3
        do a = 1, 3
          call comms_scatterv(CC_q_loc, w*counts(mpirank(comm)), CC_q(:, :, :, a, b), 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, CC_q_loc, CC_R_temp)
          call comms_reduce(CC_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, CC_R_temp, CC_R(:, :, :, a, b))
          end if
        end do
      end do

      deallocate (CC_q_loc)
      deallocate (CC_q)
      deallocate (CC_R_temp)
    end if

    call comms_bcast(CC_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_CC_R', timer)
    return

105 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.uHu', comm)
    return
106 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.uHu', comm)
    return

  end subroutine get_CC_R