pw90common_fourier_R_to_k_vec_dadb Subroutine

  subroutine pw90common_fourier_R_to_k_vec_dadb(kpt, OO_R, OO_da, OO_dadb)
    !====================================================================!
    !                                                                    !
    !! For $$OO_{ij;dx,dy,dz}$$:
    !! $$O_{ij;dx,dy,dz}(k) = \sum_R e^{+ik.R} O_{ij;dx,dy,dz}(R)$$
    !! For $$OO_{ij;dx1,dy1,dz1;dx2,dy2,dz2}$$:
    !! $$O_{ij;dx1,dy1,dz1;dx2,dy2,dz2}(k) = \sum_R e^{+ik.R} i.R_{dx2,dy2,dz2}
    !!                                       .O_{ij;dx1,dy1,dz1}(R)$$
    !                                                                    !
    !====================================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_i, twopi
    use w90_parameters, only: num_kpts, kpt_latt, num_wann, use_ws_distance
    use w90_ws_distance, only: irdist_ws, crdist_ws, wdist_ndeg, ws_translate_dist

    implicit none

    ! Arguments
    !
    real(kind=dp)                                     :: kpt(3)
    complex(kind=dp), dimension(:, :, :, :), intent(in)  :: OO_R
    complex(kind=dp), optional, dimension(:, :, :), intent(out)     :: OO_da
    complex(kind=dp), optional, dimension(:, :, :, :), intent(out)   :: OO_dadb

    integer          :: ir, i, j, ideg, a, b
    real(kind=dp)    :: rdotk
    complex(kind=dp) :: phase_fac

    if (use_ws_distance) CALL ws_translate_dist(nrpts, irvec)
    if (present(OO_da)) OO_da = cmplx_0
    if (present(OO_dadb)) OO_dadb = cmplx_0
    do ir = 1, nrpts
! [lp] Shift the WF to have the minimum distance IJ, see also ws_distance.F90
      if (use_ws_distance) then
        do j = 1, num_wann
        do i = 1, num_wann
          do ideg = 1, wdist_ndeg(i, j, ir)

            rdotk = twopi*dot_product(kpt(:), real(irdist_ws(:, ideg, i, j, ir), dp))
            phase_fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ndegen(ir)*wdist_ndeg(i, j, ir), dp)
            if (present(OO_da)) then
              OO_da(i, j, 1) = OO_da(i, j, 1) + phase_fac*OO_R(i, j, ir, 1)
              OO_da(i, j, 2) = OO_da(i, j, 2) + phase_fac*OO_R(i, j, ir, 2)
              OO_da(i, j, 3) = OO_da(i, j, 3) + phase_fac*OO_R(i, j, ir, 3)
            endif
            if (present(OO_dadb)) then
              do a = 1, 3
                do b = 1, 3
                  OO_dadb(i, j, a, b) = OO_dadb(i, j, a, b) + &
                                        cmplx_i*crdist_ws(b, ideg, i, j, ir)*phase_fac*OO_R(i, j, ir, a)
                enddo
              enddo
            endif

          enddo
        enddo
        enddo
      else
! [lp] Original code, without IJ-dependent shift:
        rdotk = twopi*dot_product(kpt(:), irvec(:, ir))
        phase_fac = cmplx(cos(rdotk), sin(rdotk), dp)/real(ndegen(ir), dp)
        if (present(OO_da)) then
          OO_da(:, :, 1) = OO_da(:, :, 1) + phase_fac*OO_R(:, :, ir, 1)
          OO_da(:, :, 2) = OO_da(:, :, 2) + phase_fac*OO_R(:, :, ir, 2)
          OO_da(:, :, 3) = OO_da(:, :, 3) + phase_fac*OO_R(:, :, ir, 3)
        endif
        if (present(OO_dadb)) then
          do a = 1, 3
            do b = 1, 3
              OO_dadb(:, :, a, b) = OO_dadb(:, :, a, b) + cmplx_i*crvec(b, ir)*phase_fac*OO_R(:, :, ir, a)
            enddo
          enddo
        endif
      endif
    enddo

  end subroutine pw90common_fourier_R_to_k_vec_dadb