SUBROUTINE mix_pot_diel (ndim, vout, vin, dr2, tr2, conv, alphamix, iter, iw)

  USE kinds,          ONLY : DP 
  USE io_global,     ONLY : stdout, ionode_id, ionode
  USE control_gw,     ONLY : nmix_gw
  USE mp_global,      ONLY : inter_pool_comm, intra_pool_comm, mpime, mp_global_end
  USE gvect,          ONLY : nrxx, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
                             nl, ngm, g, nlm
  USE constants,      ONLY : e2, fpi, RYTOEV, pi, eps8
  USE cell_base,      ONLY : alat, tpiba2, omega
  USE qpoint,         ONLY : xq
  USE mp,             ONLY : mp_barrier, mp_bcast, mp_sum
  USE freq_gw,        ONLY : fpol, fiu, nfs, nfsmax, nwcoul, wcoul
  USE gwsigma,        ONLY : ngmsco, sigma, sigma_g, nrsco, nlsco, fft6_g2r, ecutsco, ngmsig,&
                            nr1sco, nr2sco, nr3sco, ngmgrn, ngmpol
  USE units_gw,      ONLY : iuncoul, iungreen, iunsigma, lrsigma, lrcoul, lrgrn, iuwfc, lrwfc
  USE io_files,      ONLY : prefix, iunigk, prefix, tmp_dir

IMPLICIT NONE

  integer          :: ndim, iter, n_iter, n, ig, ndimtot, iw, igp
  integer*8 :: unf_recl
  complex(kind=DP) :: vout (ndim), vin (ndim), vf(ndim)
  COMPLEX(DP), ALLOCATABLE :: scrcoul_g (:,:,:)
  real(kind=DP)    :: alphamix, dr2, tr2
  logical          :: conv
  REAL(DP)         :: wwp, eps0, q0, wwq, fac, qg
  REAL(DP)         :: qg2, qg2coul
  complex(kind=DP) :: diel 
  real(kind=DP) :: DZNRM2
  complex(kind=DP) :: ZDOTC
  INTEGER :: rec0, ios
  external ZDOTC, DZNRM2
  character(len=256) :: tempfile, filename
#define DIRECT_IO_FACTOR 8 

  ! On input :
  !    ndim      dimension of arrays vout, vin
  !    vout      output potential/rho at current iteration
  ! On output :
  !    tr2       threshold for selfconsistency
  !    vin       potential/rho at previous iteration
  !    conv      true if dr2.le.tr2
  !    dr2       [(vout-vin)/ndim]^2

!Si parameters

  ALLOCATE ( scrcoul_g      (ngmpol, ngmpol, nfs)     )
  wwp    = 18.0/RYTOEV    ! plasma frequency in Ry
  eps0   = 11.4           ! static diel constant of Si
  q0     = 1.1            ! characteristic momentum of Si Resta
  fac    = 1.d0/(1.d0-1.d0/eps0)

#ifdef __PARA
if(.not.ionode) then
!OPEN coulomb file (only written to by head node).
   filename = trim(prefix)//"."//"coul1"
   tempfile = trim(tmp_dir) // trim(filename)
   unf_recl = DIRECT_IO_FACTOR * int(lrcoul, kind=kind(unf_recl))
   open(iuncoul, file = trim(adjustl(tempfile)), iostat = ios, &
   form = 'unformatted', status = 'OLD', access = 'direct', recl = unf_recl)
endif
#endif

CALL davcio(scrcoul_g, lrcoul, iuncoul, 1, -1)

!check it's actually reading scrcoul_g.
!         write(1000+mpime,'("realspace")')
!         do ig = 1, 7
!             write(1000+mpime,'(7f14.7)')real(vout(ig))
!         enddo
!         write(1000+mpime,*)
!         do ig = 1, 7
!             write(1000+mpime,'(7f14.7)')real(scrcoul_g(ig,1:7,1))
!         enddo
!MoS2 parameters
!  wwp    = 24.0/RYTOEV   ! plasma frequency in Ry
!  eps0   = 7.4           ! static diel constant of Si
!  q0     = 1.90          ! characteristic momentum of MoS2 calculated from Resta..
!Test convergence
  do n = 1, ndim
     vf (n) = vout (n) - vin (n)
  enddo

  dr2 = DZNRM2 (ndim, vf, 1)**2
  ndimtot = ndim
  dr2 = (sqrt (dr2) / ndimtot)**2
  conv = dr2.lt.tr2

!Flip (V_{out} - V_{in}) to G space.
! call cft3 (vout, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
  call cft3 (vf, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)

!do FFT properly...
  do ig = 1, ngm
     vout(ig) = vf(nl(ig))
  enddo

  do ig = 1, ngmpol
     !add one back to the diagonal since I store inveps-1
     scrcoul_g(ig,ig,iw) =  scrcoul_g(ig,ig,iw) + dcmplx(1.0d0,0.0d0)
  enddo

  vf(:) = dcmplx(0.0d0,0.0d0)

  do ig = 1, ngmpol
   do igp = 1, ngmpol
      vf(ig) = vf(ig) + scrcoul_g(ig,igp,iw)*vout(igp)
   enddo
  enddo

!scale higher elements. 
!do ig = ngmpol+1,ngm
! qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
! qg     = sqrt(tpiba2*qg2)
! wwq    = wwp * sqrt (fac * (1.d0 + (qg/eps0/q0)**2.d0))
! diel = DCMPLX(1 - wwp**2.d0/((fiu(iw))**2.d0 + wwq**2.d0), 0.d0)
! diel = DCMPLX(1 - wwp**2.d0/(wwq**2.d0), 0.d0)
! vout(ig) =  diel*vout(ig) 
!enddo

!Linea mix higher elements.
  vout(1:ngmpol) =  vin(1:ngmpol) + vf(1:ngmpol)

  do n = ngmpol+1, ngm
     vout(n) = dcmplx(0.9, 0.d0)*vin(n) + dcmplx(0.1, 0.0d0)*vout(n)
  enddo

!again do the proper fft.
  vf(:) = dcmplx(0.0d0, 0.0d0)
  do ig = 1, ngmpol 
     vf(nl(ig)) =  vout(ig)
  enddo

  call cft3 (vf, nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)
  
  vin = vf 
  vout = vin
 
  END SUBROUTINE mix_pot_diel
!-----------------------------------------------------------------------