!------------------------------------------------------
  subroutine qmatrixd(c0, bec0,ctable, gqq, qmat, detq)
!------------------------------------------------------
!
! this subroutine computes the inverse of the matrix Q
! Q_ij=<Psi_i^0|e^iG_ipol.r|Psi_j^0>
! and det Q
! respect to vectorial (serial) program I changed ngwx to ngw :-)
!    c0 input: the unperturbed wavefunctions
!    bec0 input: the coefficients <Phi_Rj|Psi_i^0>
!    ctable input: the coorespondence array
!    gqq input: the intqq(r) exp(iG_ipol*r) array
!    qmat output: the inverse q matrix
!    detq output: det Q

  use gvec
  use gvecs
  use elct
  use cnst
  use cvan
  use parm
  use parms
  use cvan
  use ions

  implicit none


  real(kind=8)      bec0(nhsa,n)
  complex(kind=8)   gqq(nhx,nhx,nas,nsp)
  complex(kind=8)   c0(ngw,nx),qmat(nx,nx),detq
  integer           ctable(ngw,2)

  complex(kind=8)   work(nx),tmp
  integer           ig,ix,jx,iv,jv,is,ia,i
  integer           ipiv(nx,nx),info,ii
  integer           ism,isa,inl,jnl

! first the local part

  do ix=1,n
     do jx=ix,n

        tmp=(0.,0.)
        do ig=1,ngw
           if(ctable(ig,1).ne.(ngw+1))then
              if(ctable(ig,1).ge.0) then
                 tmp=tmp+conjg(c0(ctable(ig,1),ix))*c0(ig,jx)
              else
                 tmp=tmp+c0(-ctable(ig,1),ix)*c0(ig,jx)
              endif
           endif
        enddo
        do ig=ng0,ngw
           if(ctable(ig,2).ne.(ngw+1)) then
              if(ctable(ig,2).lt.0) then
                 tmp=tmp+c0(-ctable(ig,2),ix)*conjg(c0(ig,jx))
              else
                 tmp=tmp+conjg(c0(ctable(ig,2),ix))*conjg(c0(ig,jx))
              endif
           endif
        enddo

        qmat(ix,jx)=tmp
        qmat(jx,ix)=tmp

     end do
  end do

  call reduce(2*nx*nx,qmat)

!  now the non local vanderbilt part

  do ix=1,n
     do jx=ix,n

        tmp=(0.,0.)
        do is=1,nvb          !loop on vanderbilt species
           do ia=1,na(is)      !loop on atoms
              do iv=1,nh(is)     !loop on projectors
                 do jv=1,nh(is)
                    inl=ish(is)+(iv-1)*na(is)+ia
                    jnl=ish(is)+(jv-1)*na(is)+ia
                    tmp=tmp+gqq(iv,jv,ia,is)*bec0(inl,ix)*bec0(jnl,jx)
                 enddo
              enddo
           enddo
        enddo

        qmat(ix,jx)=qmat(ix,jx)+tmp
        qmat(jx,ix)=qmat(ix,jx)

     enddo
  enddo



! LAPACK inverse and determinant

  call zgetrf(n,n,qmat,nx,ipiv,info)!ATTENZIONE
  detq=(1.,0.)
  do ii=1,n
     if(ii.ne.ipiv(ii,1)) detq=-detq
  enddo
  do ii=1,n
     detq = detq*qmat(ii,ii)
  enddo
  call zgetri(n,qmat,nx,ipiv,work,nx,info)

  if (detq.eq.0) then
     write(6,*) 'warning! detq = 0, replacing with 0.001'
     detq=0.001
  end if

! now qmat is symetrized

  do ix=1,n
     do jx=ix,n
        qmat(jx,ix)=0.5*(qmat(ix,jx)+qmat(jx,ix))
        qmat(ix,jx)=qmat(jx,ix)
     enddo
  enddo

  return
  end subroutine qmatrixd