!  this subroutine computes the array gqq and gqqm
!  gqq=int_dr qq(r)exp(iGr)=<Beta_r|exp(iGr)|Beta_r'>
!  gqqm=int_dr qq(r)exp(-iGr)=<Beta_r|exp(-iGr)|Beta_r'>
! ATTENZIONE ora solo cella cubica
subroutine qqberry2( gqq,gqqm, ipol)

!   gqq output: as defined above

  use gvec
  use cvan
  use core
  use cnst
  use parm
  use ncprm
  use gvecb
  use parmb
  use elct
  use ions
  use para_mod
  
  implicit none


  complex(kind=8) gqq(nhx,nhx,nas,nsp)
  complex(kind=8) gqqm(nhx,nhx,nas,nsp)
  real(kind=8) gmes
  integer ipol


! local variables

  integer ig, is, iv, jv, i, istart, il,l,ir,    & 
       &     igi,ia
  real(kind=8), allocatable:: fint(:),jl(:)
  real(kind=8), allocatable:: qradb2(:,:,:,:) 
  real(kind=8) c, ylmr, xg
  complex(kind=8) qgbs,sig
  integer ivs, jvs, ivl, jvl,   lp
  real(kind=8), allocatable:: ylm(:)
  external ylmr

  allocate( fint( mmaxx))
  allocate( jl(mmaxx))
  allocate( qradb2(nbrx,nbrx,lx,nsp))
  allocate( ylm(ngw))


  call zero(nbrx*nbrx*lx*nsp,qradb2)
     

! qui deve trovare corrispondenza ipol, g adatto

     
  do is=1,nsp
     do ia=1,nas
        do jv=1,nhx
           do iv=1,nhx
              gqq(iv,jv,ia,is)=(0.,0.)
              gqqm(iv,jv,ia,is)=(0.,0.)
           enddo
        enddo
     enddo
  enddo
  

  if(ipol.eq.1) then
     gmes=a1(1)**2+a1(2)**2+a1(3)**2
     gmes=2*pi/SQRT(gmes)
  endif
  if(ipol.eq.2) then
     gmes=a2(1)**2+a2(2)**2+a2(3)**2
     gmes=2*pi/SQRT(gmes)
  endif
  if(ipol.eq.3) then
     gmes=a3(1)**2+a3(2)**2+a3(3)**2
     gmes=2*pi/SQRT(gmes)
  endif    

  do is=1,nvb!only for Vanderbilt species 
     c=fpi                 !/omegab
! now the radial part
     do l=1,nqlc(is)                        
        xg= gmes
        call bess(xg,l,kkbeta(is),r(1,is),jl)
        do iv= 1,nbeta(is)
           do jv=iv,nbeta(is)
!     
!     note qrl(r)=r^2*q(r)
!
              do ir=1,kkbeta(is)
                 fint(ir)=qrl(ir,iv,jv,l,is)*jl(ir)
              end do
              if (ipp(is).eq.0) then
                 call herman_skillman_int     &
                      &  (kkbeta(is),cmesh(is),fint,qradb2(iv,jv,l,is))
              else
                 call simpson     &
                      &           (kkbeta(is),fint,rab(1,is),qradb2(iv,jv,l,is))                         
              endif
              qradb2(iv,jv,l,is)=  c*qradb2(iv,jv,l,is)
              qradb2(jv,iv,l,is)=  qradb2(iv,jv,l,is)
           end do
        end do
     end do
         
  enddo
  igi=-1
  do ig=1,ngw
     if(ipol.eq.1 ) then
        if(in1p(ig).eq.1 .and. in2p(ig).eq.0  .and. in3p(ig).eq. 0) igi=ig
     endif
     if(ipol.eq.2 ) then
        if(in1p(ig).eq.0 .and. in2p(ig).eq.1  .and. in3p(ig).eq. 0) igi=ig
     endif
     if(ipol.eq.3 ) then
        if(in1p(ig).eq.0 .and. in2p(ig).eq.0   .and. in3p(ig).eq. 1) igi=ig
     endif
  enddo
  if( igi.ne.-1) then

!setting array beigr

      
             
     do is=1,nvb
        do iv= 1,nh(is)
           do jv=iv,nh(is)
              ivs=indv(iv,is)
              jvs=indv(jv,is)
              ivl=indlm(iv,is)
              jvl=indlm(jv,is)
!
!     lpx = max number of allowed y_lm
!     lp  = composite lm to indentify them
!
              qgbs=(0.,0.)
              do i=1,lpx(ivl,jvl)
                 lp=lpl(ivl,jvl,i)
!
!     extraction of angular momentum l from lp:  
!
                 if (lp.eq.1) then
                    l=1         
                 else if ((lp.ge.2) .and. (lp.le.4)) then
                    l=2
                 else if ((lp.ge.5) .and. (lp.le.9)) then
                    l=3
                 else if ((lp.ge.10).and.(lp.le.16)) then
                    l=4
                 else if ((lp.ge.17).and.(lp.le.25)) then
                    l=5
                 else if (lp.ge.26) then 
                    call error(' qvanb ',' lp.ge.26 ',lp)
                 endif
!
!       sig= (-i)^l
!
                 sig=(0.,-1.)**(l-1)
                  
!                  call ylmr2b(lp,4,gx,ylm) ATTENZIONE
!                  qgbs=qgbs+sig*ylm(igi)*qradb2(2,ivs,jvs,l,is)
                 sig=sig*ap(lp,ivl,jvl)
                 qgbs=qgbs+sig*ylmr(lp,igi)*qradb2(ivs,jvs,l,is)
                
              end do
   

              
              do ia=1,na(is)
!                  gqq(iv,jv,ia,is)=qgbs*eigr(igi,ia,is)!ATTENZIONE era cosi'
!                  gqq(jv,iv,ia,is)=qgbs*eigr(igi,ia,is)
!                  gqqm(iv,jv,ia,is)=conjg(gqq(iv,jv,ia,is))
!                  gqqm(jv,iv,ia,is)=conjg(gqq(iv,jv,ia,is))
                     
                 gqqm(iv,jv,ia,is)=qgbs
                 gqqm(jv,iv,ia,is)=qgbs
                 gqq(iv,jv,ia,is)=conjg(gqqm(iv,jv,ia,is))
                 gqq(jv,iv,ia,is)=conjg(gqqm(iv,jv,ia,is))
              end do
           end do
        enddo
     enddo
  endif

#ifdef PARA
  call reduce(2*nhx*nhx*nas*nsp, gqq)
  call reduce(2*nhx*nhx*nas*nsp, gqqm)
#endif


  deallocate( fint)
  deallocate( jl)
  deallocate(qradb2)
  deallocate(ylm)
  
  return
end subroutine qqberry2






! this subroutine updates gqq and gqqm to the 
! (new) atomic position


subroutine qqupdate(eigr, gqqm0, gqq, gqqm, ipol)

!   gqq output: as defined above

  use cvan
  use ions
  use gvec
  use elct

  implicit none

 
  complex(kind=8) eigr(ngw,nas,nsp)
  complex(kind=8) gqq(nhx,nhx,nas,nsp)
  complex(kind=8) gqqm(nhx,nhx,nas,nsp)
  complex(kind=8) gqqm0(nhx,nhx,nas,nsp)

  integer ipol
  
  integer igi,ig,is,iv,jv,ia


  do is=1,nsp
     do ia=1,nas
        do jv=1,nhx
           do iv=1,nhx
              gqq(iv,jv,ia,is)=(0.,0.)
              gqqm(iv,jv,ia,is)=(0.,0.)
           enddo
        enddo
     enddo
  enddo

  igi=-1
  do ig=1,ngw
     if(ipol.eq.1 ) then
        if(in1p(ig).eq.1 .and. in2p(ig).eq.0  .and. in3p(ig).eq. 0) igi=ig
     endif
     if(ipol.eq.2 ) then
        if(in1p(ig).eq.0 .and. in2p(ig).eq.1  .and. in3p(ig).eq. 0) igi=ig
     endif
     if(ipol.eq.3 ) then
        if(in1p(ig).eq.0 .and. in2p(ig).eq.0   .and. in3p(ig).eq. 1) igi=ig
     endif
  enddo
  if( igi.ne.-1) then

  
     do is=1,nvb
        do ia=1,na(is)
           do iv= 1,nh(is)
              do jv=iv,nh(is)
                 gqqm(iv,jv,ia,is)= gqqm0(iv,jv,ia,is)*eigr(igi,ia,is)
                 gqqm(jv,iv,ia,is)= gqqm0(iv,jv,ia,is)*eigr(igi,ia,is)
                 gqq(iv,jv,ia,is)=conjg(gqqm(iv,jv,ia,is))
                 gqq(jv,iv,ia,is)=conjg(gqqm(iv,jv,ia,is))
              enddo
           enddo
        enddo
     enddo
  endif
#ifdef PARA
  call reduce(2*nhx*nhx*nas*nsp, gqq)
  call reduce(2*nhx*nhx*nas*nsp, gqqm)
#endif
  return
end subroutine qqupdate