pw90common_fourier_R_to_k_new_second_d Subroutine

  subroutine pw90common_fourier_R_to_k_new_second_d(kpt, OO_R, OO, OO_da, OO_dadb)
    !=======================================================!
    !                                                       !
    !! For OO:
    !! $$O_{ij}(k) = \sum_R e^{+ik.R}.O_{ij}(R)$$
    !! For $$OO_{dx,dy,dz}$$:
    !! $$\sum_R [i.R_{dx,dy,dz}.e^{+ik.R}.O_{ij}(R)]$$
    !! where R_{x,y,z} are the Cartesian components of R
    !! For $$OO_{dx1,dy1,dz1;dx2,dy2,dz2}$$:
    !! $$-\sum_R [R_{dx1,dy1,dz1}.R_{dx2,dy2,dz2}.e^{+ik.R}.O_{ij}(R)]$$
    !! where R_{xi,yi,zi} are the Cartesian components of R
    !                                                       !
    !=======================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_i, twopi
    use w90_parameters, only: timing_level, 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
    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)) OO = cmplx_0
    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)) OO(i, j) = OO(i, j) + phase_fac*OO_R(i, j, ir)
            if (present(OO_da)) then
              do a = 1, 3
                OO_da(i, j, a) = OO_da(i, j, a) + cmplx_i*crdist_ws(a, ideg, i, j, ir)* &
                                 phase_fac*OO_R(i, j, ir)
              enddo
            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) - crdist_ws(a, ideg, i, j, ir)* &
                                        crdist_ws(b, ideg, i, j, ir)*phase_fac*OO_R(i, j, ir)
                enddo
              enddo
            end if

          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)) OO(:, :) = OO(:, :) + phase_fac*OO_R(:, :, ir)
        if (present(OO_da)) then
          do a = 1, 3
            OO_da(:, :, a) = OO_da(:, :, a) + cmplx_i*crvec(a, ir)*phase_fac*OO_R(:, :, ir)
          enddo
        endif
        if (present(OO_dadb)) then
          do a = 1, 3
            do b = 1, 3
              OO_dadb(:, :, a, b) = OO_dadb(:, :, a, b) - &
                                    crvec(a, ir)*crvec(b, ir)*phase_fac*OO_R(:, :, ir)
            enddo
          enddo
        end if
      endif
    enddo

  end subroutine pw90common_fourier_R_to_k_new_second_d