!
  !-----------------------------------------------------------------------
  subroutine readmat_shuffle2 ( iq_irr, nqc_irr, nq, iq_first, sxq, imq)
  !-----------------------------------------------------------------------
  !
  ! read dynamical matrix for the q points  
  ! iq_first, iq_first+1, ... iq_first+nq-1
  !
  !-----------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  use io_files, only : prefix, tmp_dir, dvscf_dir
  USE cell_base, ONLY : ibrav, celldm, omega
  USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
  USE el_phon, ONLY : dynq
  !
  implicit none
  !
  ! Input
  !
  integer :: iq_irr, nqc_irr, nq, iq_first, imq
  !  the index of the irreducible q point 
  !  the number of irreducible qpoints 
  !  the number of q points in the star of q
  !  the index of the first qpoint to be read in the uniform q-grid
  !  flag which tells whether we have to consider the -q vectors
  real(kind=DP) :: sxq(3,48)
  !  the q vectors in the star
  !
  ! output
  !
  ! dynq (nmode,nmodes,nqc) (in el_phon.mod)
  !
  ! local 
  !
  real(kind=DP) :: eps, celldm_ (6), amass_, tau_ (3), q(3), &
       dynr (2, 3, nat, 3, nat), massfac, sumr(2)
  integer :: iudyn, ntyp_, nat_, ibrav_, ityp_, ios, iq, &
       nt, na, nb, naa, nbb, nu, mu, i, j, ipol, jpol
  character(len=80) :: line
  character(len=3)  :: atm, filelab
  character (len=256) :: tempfile
  !
  eps = 1.d-5
  !
  !  the call to set_ndnmbr is just a trick to get quickly
  !  a file label by exploiting an existing subroutine
  !  (if you look at the sub you will find that the original
  !  purpose was for pools and nodes)
  !
  call set_ndnmbr ( 0, iq_irr, 1, nqc_irr, filelab)
  tempfile = trim(dvscf_dir) // trim(prefix) // '.dyn_q' // filelab
  !
  iudyn = 81
  open (unit = iudyn, file = tempfile, status = 'unknown', err = 100, iostat = ios)
100  call errore ('openfilq', 'opening file'//tempfile, abs (ios) )
  rewind (iudyn)
  !
  !  read header and run some checks
  !
  read (iudyn, '(a)') line
  read (iudyn, '(a)') line
  read (iudyn, * ) ntyp_, nat_, ibrav_, celldm_
  if (ntyp.ne.ntyp_.or.nat.ne.nat_.or.ibrav_.ne.ibrav.or.abs ( &
     celldm_ (1) - celldm (1) ) .gt.1.0d-5) call errore ('readmat', &
     'inconsistent data', 1)
  do nt = 1, ntyp
     read (iudyn, * ) i, atm, amass_
     if (nt.ne.i.or.abs (amass_ - amass (nt) ) .gt.1.0d-5) call errore ( &
          'readmat', 'inconsistent data', 1 + nt)
  enddo
  do na = 1, nat
     read (iudyn, * ) i, ityp_, tau_
     if (na.ne.i.or.ityp_.ne.ityp (na) ) call errore ('readmat', &
        'inconsistent data', 10 + na)
  enddo
  !
  !  read dyn mat for all q in the star
  !
  do iq = 1, nq
    !
    read (iudyn, '(///a)') line
    read (line (11:80), * ) (q (i), i = 1, 3)
    !
    ! if imq=0 we have to compare with -sxq
    !
    if (imq.ne.0) then
      if ( abs(q(1)-sxq(1,iq)).gt.eps .or. &
           abs(q(2)-sxq(2,iq)).gt.eps .or. &
           abs(q(3)-sxq(3,iq)).gt.eps ) call errore &
         ('readmat', 'wrong qpoint', 0)
    else
      if ( abs(q(1)+sxq(1,iq)).gt.eps .or. &
           abs(q(2)+sxq(2,iq)).gt.eps .or. &
           abs(q(3)+sxq(3,iq)).gt.eps ) call errore &
         ('readmat', 'wrong qpoint', 0)
    endif
    read (iudyn, '(a)') line
    !
    !
    do na = 1, nat
       do nb = 1, nat
          read (iudyn, * ) naa, nbb
          if (na.ne.naa.or.nb.ne.nbb) call errore &
               ('readmat', 'error reading file', nb)
          read (iudyn, * ) &
              ((dynr (1,i,na,j,nb), dynr (2,i,na,j,nb), j = 1, 3), i = 1, 3)
       enddo
    enddo
    !
    if ( abs(q(1)).lt.eps .and. abs(q(2)).lt.eps .and. abs(q(3)).lt.eps ) then
      !
      ! in the case q=0 we impose the acoustic sum rule
      ! [Gonze and Lee, PRB 55, 10361 (1998), Eq. (45)]
      !
      write(6,'(8x,a)') 'Imposing acoustic sum rule on the dynamical matrix'
      !
      do na = 1,nat
       do ipol = 1,3
        do jpol = ipol,3
          !
          sumr(1) = sum ( dynr (1,ipol,na,jpol,:) )
          sumr(2) = sum ( dynr (2,ipol,na,jpol,:) )
          !
          dynr (:,ipol,na,jpol,na) = dynr (:,ipol,na,jpol,na) - sumr
          !
        end do
       end do
      end do
      !
    endif
    !
    !  fill the two-indices dynamical matrix in cartesian coordinates
    !  the proper index in the complete list is iq_first+iq-1
    !
    do na = 1,nat
     do nb = 1,nat
        do ipol = 1,3
         do jpol = 1,3
           !
           mu = (na-1)*3+ipol
           nu = (nb-1)*3+jpol
           dynq ( mu, nu, iq_first+iq-1) = &
             dcmplx ( dynr (1,ipol,na,jpol,nb), dynr (2,ipol,na,jpol,nb) )
           !
         end do
        end do
     end do
    end do
    !
  enddo
  !
  close(iudyn)
  !
  return
  end subroutine readmat_shuffle2
  !