!
! Copyright (C) 2001-2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
subroutine dvqpsi_us (igpert, ik, uact, addnlcc)
  !----------------------------------------------------------------------
  !
  ! This routine calculates dV_bare/dtau * psi for one perturbation
  ! with a given q. The displacements are described by a vector u.
  ! The result is stored in dvpsi. The routine is called for each k point
  ! and for each pattern u. It computes simultaneously all the bands.
  ! It implements Eq. B29 of PRB 64, 235118 (2001). The contribution
  ! of the local pseudopotential is calculated here, that of the nonlocal
  ! pseudopotential in dvqpsi_us_only.
  !
  ! HL- FORMERLY this routine does what it says up there. I've commented out the stuff I think pertains
  ! to initializing a phonon perturbation. Then stuck the igpert counter in so that I can switch on 
  ! one fourier component of the perturbing potential, then FFT it to real space, and continue on the way.
  !

  USE kinds, only : DP
  USE ions_base, ONLY : nat, ityp
  USE cell_base, ONLY : tpiba
  USE gvect,     ONLY : nrxx, eigts1, eigts2, eigts3, ig1,ig2,ig3, g, nl, &
                        ngm, nr1,nr2,nr3,nrx1,nrx2,nrx3
  USE gsmooth,   ONLY : nrxxs, ngms, doublegrid, nls, &
                        nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
  USE lsda_mod,  ONLY : lsda, isk
  USE noncollin_module, ONLY : npol
  use uspp_param,ONLY : upf
  USE wvfct,     ONLY : nbnd, npw, npwx, igk
  USE wavefunctions_module,  ONLY: evc
  USE nlcc_gw,    ONLY : nlcc_any, drc
  USE eqv,        ONLY : dvpsi, dmuxc, vlocq
  USE qpoint,     ONLY : npwq, igkq, xq, eigqts, ikks

  implicit none
  !
  !   The dummy variables
  !

  integer :: ik, igpert
  ! input: the k point
  ! counter for G-vector for plane wave perturbation  

  complex(DP) :: uact (3 * nat)
  ! input: the pattern of displacements

  logical :: addnlcc
  !
  !   And the local variables
  !

  integer :: na, mu, ikk, ig, nt, ibnd, ir, is, ip
  ! counter on atoms
  ! counter on modes
  ! the point k
  ! counter on G vectors
  ! the type of atom
  ! counter on bands
  ! counter on real mesh

  complex(DP) :: gtau, gu, fact, u1, u2, u3, gu0
  complex(DP) , allocatable, target :: aux (:)
  complex(DP) , allocatable :: aux1 (:), aux2 (:)
  complex(DP) , pointer :: auxs (:)
  ! work space

  call start_clock ('dvqpsi_us')
  if (nlcc_any.and.addnlcc) then
     allocate (aux( nrxx))    
     if (doublegrid) then
        allocate (auxs( nrxxs))    
     else
        auxs => aux
     endif
  endif

  allocate (aux1( nrxxs))    
  allocate (aux2( nrxxs))    
  !
  !    We start by computing the contribution of the local potential.
  !    The computation of the derivative of the local potential is done in
  !    reciprocal space while the product with the wavefunction is done in
  !    real space
  !
  ikk = ikks(ik)
  dvpsi(:,:) = (0.d0, 0.d0)
  aux1(:) = (0.d0, 0.d0)

#if 0
  do na = 1, nat
     fact = tpiba * (0.d0, -1.d0) * eigqts (na)
     mu = 3 * (na - 1)
     if (abs (uact (mu + 1) ) + abs (uact (mu + 2) ) + abs (uact (mu + &
          3) ) .gt.1.0d-12) then
        nt = ityp (na)
        u1 = uact (mu + 1)
        u2 = uact (mu + 2)
        u3 = uact (mu + 3)
        gu0 = xq (1) * u1 + xq (2) * u2 + xq (3) * u3
        do ig = 1, ngms
           gtau = eigts1 (ig1 (ig), na) * eigts2 (ig2 (ig), na) * eigts3 ( &
                ig3 (ig), na)
           gu = gu0 + g (1, ig) * u1 + g (2, ig) * u2 + g (3, ig) * u3
           aux1 (nls (ig) ) = aux1 (nls (ig) ) + vlocq (ig, nt) * gu * &
                fact * gtau
        enddo
     endif
  enddo
#endif

!one plane wave HL
!aux1 (nls (igpert) ) = (1.d0, 0.d0)
!write (6, '("igpert ", i4)' )igpert

aux1 ( nl (igpert) ) = (1.d0, 0.d0)

!do ig = 1 , 10
!write(6,'("vlocq =  ", 3f7.3)') vlocq(ig, 1)
!enddo


!
! add NLCC when present
!
   if (nlcc_any.and.addnlcc) then
      aux(:) = (0.d0, 0.d0)
      do na = 1,nat
         fact = tpiba*(0.d0,-1.d0)*eigqts(na)
         mu = 3*(na-1)
         if (abs(uact(mu+1))+abs(uact(mu+2))  &
                         +abs(uact(mu+3)).gt.1.0d-12) then
            nt=ityp(na)
            u1 = uact(mu+1)
            u2 = uact(mu+2)
            u3 = uact(mu+3)
            gu0 = xq(1)*u1 +xq(2)*u2+xq(3)*u3
            if (upf(nt)%nlcc) then
               do ig = 1,ngm
                  gtau = eigts1(ig1(ig),na)*   &
                         eigts2(ig2(ig),na)*   &
                         eigts3(ig3(ig),na)
                  gu = gu0+g(1,ig)*u1+g(2,ig)*u2+g(3,ig)*u3
                  aux(nl(ig))=aux(nl(ig))+drc(ig,nt)*gu*fact*gtau
               enddo
            endif
         endif
      enddo
      call cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
      if (.not.lsda) then
         do ir=1,nrxx
            aux(ir) = aux(ir) * dmuxc(ir,1,1)
         end do
      else
         is=isk(ikk)
         do ir=1,nrxx
            aux(ir) = aux(ir) * 0.5d0 *  &
                 (dmuxc(ir,is,1)+dmuxc(ir,is,2))
         enddo
      endif

     call cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)

      if (doublegrid) then
         auxs(:) = (0.d0, 0.d0)
         do ig=1,ngms
            auxs(nls(ig)) = aux(nl(ig))
         enddo
      endif
      aux1(:) = aux1(:) + auxs(:)
   endif

  !
  !PH- Now we compute dV_loc/dtau in real space
  !HL- FFT dvbare into R space and compute Delta V_q(r')*Psi_nk(r') 
  !

  call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 1)

  do ibnd = 1, nbnd
     do ip=1,npol
        aux2(:) = (0.d0, 0.d0)
        if (ip==1) then
           do ig = 1, npw
              aux2 (nls (igk (ig) ) ) = evc (ig, ibnd)
           enddo
        else
           do ig = 1, npw
              aux2 (nls (igk (ig) ) ) = evc (ig+npwx, ibnd)
           enddo
        end if
        !
        !  This wavefunction is computed in real space
        !

        call cft3s (aux2, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
        do ir = 1, nrxxs
           aux2 (ir) = aux2 (ir) * aux1 (ir)
        enddo

        !
        ! and finally dV_loc/dtau * psi is transformed in reciprocal space
        !

        call cft3s (aux2, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)

        if (ip==1) then
           do ig = 1, npwq
              dvpsi (ig, ibnd) = aux2 (nls (igkq (ig) ) )
           enddo
        else
           do ig = 1, npwq
              dvpsi (ig+npwx, ibnd) = aux2 (nls (igkq (ig) ) )
           enddo
        end if
     enddo
  enddo
!
!HL FFT how it works? 
!This should recover the perturbing potential in G space. 
!         write (6, '("")' )
!         do ig = 1, 3
!         write (6, '("dvbare(nl(igpert)) ", 3f7.4, 3f7.3)' )aux1(nl(ig))
!         end do
!         call cft3 (aux1, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
!         call cft3 (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -1)
!         write(6, '("")')
!         do ig = 1, 3
!         write (6, '("dvbare(nl(igpert)) ", 3f7.4, 3f7.3)' )aux1(nl(ig))
!         write (6, '("dvbare(nl(igpert)) ", 3f7.4, 3f7.3)' )aux1(ig)
!         end do
!HL This does recover the correct array of fourier components. I think drho should transform in exactly the same way.
!

  deallocate (aux2)
  deallocate (aux1)
  if (nlcc_any.and.addnlcc) then
     deallocate (aux)
     if (doublegrid) deallocate (auxs)
  endif

!  We add the contribution of the nonlocal potential in the US form
!  First a term similar to the KB case.
!  Then a term due to the change of the D coefficients. 
!  HL I don't think we need to add anything else to the right hand side.
!  No change in Q, Beta, or S of the self-consistent potential. Just the shift due to 
!  electrons. However there may be a residual term due to the NLCC when form the pseudopotential which
! needs to be added... Solution: Don't use NLCC for the time being...

!HL call dvqpsi_us_only (ik, uact)

  call stop_clock ('dvqpsi_us')
  return
end subroutine dvqpsi_us