pw90common_fourier_R_to_k_new_second_d_TB_conv Subroutine

  subroutine pw90common_fourier_R_to_k_new_second_d_TB_conv(kpt, OO_R, oo_a_R, num_wann, &
                                                            ws_region, wannier_data, real_lattice, &
                                                            mp_grid, ws_distance, wigner_seitz, &
                                                            error, comm, OO, OO_da, OO_dadb)
    !================================================!
    ! modified version of pw90common_fourier_R_to_k_new_second_d, includes wannier centres in
    ! the exponential inside the sum (so called TB convention)
    !
    !! For OO:
    !! $$O_{ij}(k) = \sum_R e^{+ik.(R+tau_ij)}.O_{ij}(R)$$
    !! For $$OO_{dx,dy,dz}$$:
    !! $$\sum_R [i.(R+tau_ij)_{dx,dy,dz}.e^{+ik.(R+tau_ij)}.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+tau_ij)_{dx1,dy1,dz1}.(R+tau_ij)_{dx2,dy2,dz2}.e^{+ik.(R+tau_ij)}.O_{ij}(R)]$$
    !! where {xi,yi,zi} denote the Cartesian components and
    !  tau_ij = tau_j - tau_i, being tau_i=<0i|r|0i> the individual wannier centres
    !================================================!

    use w90_constants, only: dp, cmplx_0, cmplx_i, twopi
    use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type
    use w90_ws_distance, only: ws_translate_dist
    use w90_utility, only: utility_cart_to_frac, utility_inverse_mat
    use w90_postw90_types, only: wigner_seitz_type
    use w90_comms, only: w90_comm_type

    implicit none

    ! arguments
    type(ws_region_type), intent(in) :: ws_region
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(in) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_wann

    real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3)

    complex(kind=dp), intent(in) :: oo_a_R(:, :, :, :)
    complex(kind=dp), intent(in) :: OO_R(:, :, :)
    complex(kind=dp), optional, intent(out) :: OO(:, :)
    complex(kind=dp), optional, intent(out) :: OO_da(:, :, :)
    complex(kind=dp), optional, intent(out) :: OO_dadb(:, :, :, :)

    ! local variables
    real(kind=dp) :: inv_lattice(3, 3)
    integer :: ir, i, j, ideg, a, b
    real(kind=dp) :: rdotk
    real(kind=dp) :: wannier_centres_frac(3, num_wann)
    real(kind=dp) :: r_sum(3)
    complex(kind=dp) :: phase_fac

    r_sum = 0.d0

    ! rotate wannier centres from cartesian to fractional coordinates
    wannier_centres_frac(:, :) = 0.d0
    call utility_inverse_mat(real_lattice, inv_lattice)
    do ir = 1, num_wann
      call utility_cart_to_frac(wigner_seitz%wannier_centres_from_AA_R(:, ir), &
                                wannier_centres_frac(:, ir), inv_lattice)
    end do

    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, wigner_seitz%nrpts_pw90
! [lp] Original code, without IJ-dependent shift:
      do j = 1, num_wann
        do i = 1, num_wann
          r_sum(:) = real(wigner_seitz%irvec_pw90(:, ir)) &
                     + wannier_centres_frac(:, j) - wannier_centres_frac(:, i)
          rdotk = twopi*dot_product(kpt(:), r_sum(:))
          phase_fac = cmplx(cos(rdotk), sin(rdotk), 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* &
                               (wigner_seitz%crvec_pw90(a, ir) + wigner_seitz%wannier_centres_from_AA_R(a, j) - &
                                wigner_seitz%wannier_centres_from_AA_R(a, i))*phase_fac*OO_R(i, j, ir)
            end do
          end if
          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) - &
                  (wigner_seitz%crvec_pw90(a, ir) + wigner_seitz%wannier_centres_from_AA_R(a, j) - &
                   wigner_seitz%wannier_centres_from_AA_R(a, i))*(wigner_seitz%crvec_pw90(b, ir) + &
                                                                  wigner_seitz%wannier_centres_from_AA_R(b, j) - &
                                                                  wigner_seitz%wannier_centres_from_AA_R(b, i))* &
                  phase_fac*OO_R(i, j, ir)
              end do
            end do
          end if
        end do
      end do
    end do

  end subroutine pw90common_fourier_R_to_k_new_second_d_TB_conv