!-------------------------------------------------------------------
  subroutine bforceion(fion,tfor,ipol,qmatinv,bec0,becdr,gqq,evalue)
!-------------------------------------------------------------------
!
! this subroutine compute the part of force for the ions due to
! electronic berry phase( see internal notes)
! it needs becdr
!
! fion       : input, forces on ions
! tfor       : input, if true it computes force
! a1,a2,a3   : input, direct lattice vectors
! ipol       : input, electric field polarization
! qmatinv    : input, inverse of Q matrix: Q_i,j=<Psi_i|exp(iG*r)|Psi_j>
! bec0       : input, factors <beta_iR|Psi_j>
! becdr      : input, factors d<beta_iR|Psi_j>/dR
! gqq        : input, Int_e exp(iG*r)*q_ijR(r)
! evalue     : input, scale of electric field

  use ions
  use cvan
  use elct
  use parm
  use cnst

  implicit none

  real(kind=8)    evalue
  complex(kind=8) qmatinv(nx,nx),gqq(nhx,nhx,nas,nsp)
  real(kind=8)    bec0(nhsa,n),becdr(nhsa,n,3)
  real(kind=8)    fion(3,nax,nsx)
  integer         ipol
  logical         tfor

  complex(kind=8) ci,tmp0,tmp1,tmp2(3),tmp3(3),gqqtmp
  complex(kind=8) utmp0_inl(n),utmp0_jnl(n),alpha0,beta0
  complex(kind=8) bec0_inl(n),bec0_jnl(n),becdr_inl(n),becdr_jnl(n)
  real(kind=8)    gmes
  integer         ix,iv,jv,ia,is,k,i,j,isa,ilm,jlm,inl,jnl,ism
  complex(kind=8) ZDOTU

  if(.not.tfor) return
  ci=(0.,1.)

  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

  alpha0=(1.,0.)
  beta0 =(0.,0.)

  do is=1,nvb
     do ia=1,na(is)

        tmp1=(0.,0.)
        do ix=1,3
           tmp3(ix)=(0.,0.)
        end do

        do iv=1,nh(is)

           inl=ish(is)+(iv-1)*na(is)+ia

           bec0_inl=bec0(inl,:)

           call ZGEMV('N',n,n,alpha0,qmatinv,n,bec0_inl,1,beta0,utmp0_inl,1)

           do jv=1,nh(is)

              jnl=ish(is)+(jv-1)*na(is)+ia

              bec0_jnl=bec0(jnl,:)

              call ZGEMV('N',n,n,alpha0,qmatinv,n,bec0_jnl,1,beta0,utmp0_jnl,1)

              tmp0=ZDOTU(n,bec0_jnl,1,utmp0_inl,1)

              tmp1=tmp1+tmp0*ci*gmes*gqq(iv,jv,ia,is)

              do ix=1,3

                 becdr_inl=becdr(inl,:,ix)
                 becdr_jnl=becdr(jnl,:,ix)

                 tmp2(ix)=(0.,0.)

                 tmp0=ZDOTU(n,becdr_inl,1,utmp0_jnl,1)

                 tmp2(ix)=tmp2(ix)+tmp0

                 tmp0=ZDOTU(n,becdr_jnl,1,utmp0_inl,1)

                 tmp2(ix)=tmp2(ix)+tmp0

              end do

              tmp3=tmp3+tmp2*gqq(iv,jv,ia,is)

           end do
        end do

        fion(ipol,ia,is)=fion(ipol,ia,is)-2.*evalue*aimag(tmp1)/gmes
        fion(:,ia,is)=fion(:,ia,is)-2.*evalue*aimag(tmp3(:))/gmes

     end do
  end do

  return
  end subroutine bforceion