!
! Copyright (C) 2001-2015 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 syme (dvtosym)
  !-----------------------------------------------------------------------
  !
  USE kinds,            only : DP
  USE constants,        ONLY : tpi
  USE fft_base,         ONLY : dfftp
  USE cell_base,        ONLY : at
  USE symm_base,        ONLY : s, ftau
  USE noncollin_module, ONLY : nspin_lsda, nspin_mag
  USE modes,            ONLY : minus_q, nsymq, irotmq, gi, gimq
  !  
  IMPLICIT NONE
  !
  COMPLEX(DP) :: dvtosym(dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, nspin_mag)
  ! the charge density response to be symmetrized
  INTEGER :: is, ri, rj, rk, i, j, k, ipol, isym, irot
  ! counters
  REAL(DP) :: gf(3), n(3)
  ! temp variables
  COMPLEX(DP), ALLOCATABLE :: dvsym(:,:,:)
  ! the symmetrized charge density response
  COMPLEX(DP) ::  aux2, term(3, 48), phase(48)
  ! auxiliary space
  ! the multiplication factor
  ! the phase factor
  !
  REAL(DP) :: sgi(3,48)
  IF (nsymq==1) RETURN
  !
  CALL start_clock ('symdvscf')
  !
  ALLOCATE(dvsym(dfftp%nr1x, dfftp%nr2x, dfftp%nr3x))
  !
  n(1) = tpi / DBLE (dfftp%nr1)
  n(2) = tpi / DBLE (dfftp%nr2)
  n(3) = tpi / DBLE (dfftp%nr3)
!!!!!!!!!!!!!!PHASE for Sq = q+G!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  DO isym = 1, nsymq
     !
     gf(:) =  gi(1,isym) * at(1, :) * n(:) + &
              gi(2,isym) * at(2, :) * n(:) + &
              gi(3,isym) * at(3, :) * n(:)
     !
     term(:,isym) = CMPLX(cos(gf(:)), sin(gf(:)), kind=DP)
     !
  ENDDO
!!!!!!!!!!!!!!PHASE for Sq = q+G!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!PHASE for S^{-1}(q+G) = q+S^{-1}G!!!!!!!!!!!!
  sgi(:,:) = 0.0d0
  DO isym = 1, nsymq
     gf(:) =  sgi(1,isym) * at(1, :) * n(:) + &
              sgi(2,isym) * at(2, :) * n(:) + &
              sgi(3,isym) * at(3, :) * n(:)
     sgterm(:,isym) = CMPLX(cos(gf(:)), sin(gf(:)), kind=DP)
  ENDDO
!!!!!!!!!!!!!!PHASE for S^{-1}(q+G) = q+S^{-1}G!!!!!!!!!!!!
  DO is = 1, nspin_lsda
     !
     dvsym(:,:,:) = (0.d0, 0.d0)
     !
     DO isym = 1, nsymq
        phase(isym) = (1.d0, 0.d0)
     ENDDO
     !
     DO k = 1, dfftp%nr3
        DO j = 1, dfftp%nr2
           DO i = 1, dfftp%nr1
              !
              ! Loop on the symmetry operations of the small group of q.
              !
              DO isym = 1, nsymq
                 !
                 ! Rotation and fractional translation: S^-1 * r - ftau
                 !
                 CALL ruotaijk (s(1,1,isym), ftau(1,isym), i, j, k, &
                 dfftp%nr1, dfftp%nr2, dfftp%nr3, ri, rj, rk)
                 !
                 ! Calculate drho(S^-1 * r - ftau) * exp(i G*r)
                 !
                 dvsym(i,j,k) = dvsym(i,j,k) + dvtosym(ri,rj,rk,is) * phase(isym)
                 !
              ENDDO
              !
              DO isym = 1, nsymq
                 phase (isym) = phase (isym) * term (1, isym)*sgterm(1,isym)
              ENDDO
              !
           ENDDO
           !
           DO isym = 1, nsymq
              phase (isym) = phase (isym) * term (2, isym)*sgterm(1,isym)
           ENDDO
           !
        ENDDO
        !
        DO isym = 1, nsymq
           phase (isym) = phase (isym) * term (3, isym)*sgterm(1,isym)
        ENDDO
        !
     ENDDO
     !Normalize
     dvtosym(:,:,:,is) = dvsym(:,:,:) / DBLE (nsymq)
  ENDDO
  !
  DEALLOCATE(dvsym)
  !
  CALL stop_clock ('symdvscf')
  !
  RETURN
  !
END SUBROUTINE syme