!-----------------------------------------------------------------------
  ! Copyright (C) 2010-2015 Henry Lambert, Feliciano Giustino
  ! 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 mod_diel(scrcoul, xq_ibk, w, screening)

  USE kinds,       ONLY : DP
  USE constants,   ONLY : e2, fpi, RYTOEV, pi, eps8
  USE gvect,       ONLY : ngm, g
  USE gwsigma,     ONLY : sigma_c_st
  USE cell_base,   ONLY : tpiba2, tpiba, omega, alat, at
  USE mp_world,      ONLY : nproc, mpime

IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: xq_ibk(3)
  REAL(DP) :: w, inveps
  REAL(DP)    :: wwp, eps0, q0, wwq, fac, z, xq(3)
  REAL(DP)    :: wwp1, wwp2, qxy, meff, invepsqg
  REAL(DP)    :: rhoi, rhoj , wwqi, wwji
  REAL(DP)    :: wwq1, wwq2
  REAL(DP)    :: A, B, wwmi, wwpl
  REAL(DP)    :: epsq, polq
  REAL(DP)    :: qg, qg2, kf2, kf, alpha
  REAL(DP)    :: expqxy, d, x 
  INTEGER     :: screening, ig
  COMPLEX(DP) :: scrcoul (sigma_c_st%ngmt, sigma_c_st%ngmt)
  LOGICAL :: testnan

!Screening selects from one of the following model dielectric functions.
!in all cases the dielectric function is diagonal in G space so we include no
!local field effects. However the model is parameterized to include the possibility
!of two modes occurring from coupling between two electron gases of different density.

!The 'magic dielectric function' Inkson 1972
!Silicon parameters...
!check resta for parameters Phys. Rev. B 16, 2717 2722 (1977)
!wwp    = 16.0/RYTOEV  ! plasma frequency in Ry
!eps0   = 11.4         ! static diel constant of Si
!q0     = 1.1          ! characteristic momentum of Si, a.u. from Resta
!LiCL parameters
!wwp    = 17.0/RYTOEV   ! plasma frequency in Ry
!eps0   = 11.04         ! static diel constant of 
!q0     = 1.2           ! characteristic momentum of Si, a.u. from Resta
!MoS2 there are two well defined excitation for parallel and perpendicular.
!this is an anisotropic material though!
!should have (at least) 2 different plasmons!
     wwp    = 7.0/RYTOEV ! plasma frequency in Ry
     eps0   = 1.01        ! static diel constant of MoS2
     q0     = 1.90        ! characteristic momentum of MoS2 calculated from Resta paper..
!!!!!!
     fac    = 1.d0/(1.d0-1.d0/eps0)
!effective mass annd density of electrons in layer i,j.
     meff   = 1.0d0
!or should this be the fourier co-efficient of the density at the particular G vector?
     rhoi   = 0.05d0
     rhoj   = 0.05d0
!z is interlayer distance.
     z      = 8.5d0
!3d: wavevector
!    kf     = (3*pi**2*rhoi)**(1/3)
!2d: wavevector
     kf     = (2*pi*rhoi)**(1/2)

 if(screening.eq.1) then
!Standard 3D-bulk ppm:
!drhoscfs (nl(ig), 1)  = 1.d0 - wwp**2.d0/((fiu(iw) + eta)**2.d0 + wwq**2.d0)
!    dvscf ( nl(ig) ) = 1.d0 - wwp**2.d0/((wim_ryd(iw)+eta)**2.d0+wwq**2.d0)
!(W-v) = (inveps(w) - delta) v
!Wave vector:
    DO ig = 1, sigma_c_st%ngmt
       qg2    = (g(1,ig)+xq_ibk(1))**2 + (g(2,ig)+xq_ibk(2))**2 + (g(3,ig)+xq_ibk(3))**2
       qg     = sqrt(tpiba2*qg2)
       fac    = 1.d0/(1.d0-(1.d0/eps0))
       wwq    = wwp * sqrt ( fac * (1.d0 + (qg/eps0/q0)**2.d0 ) )
       scrcoul(ig,ig) =  - wwp**2.d0/(w**2.d0 + wwq**2.d0)
    ENDDO
 else if(screening.eq.2) then
 !Effective bi-layer plasmon model
 !Inkson and White semicond. sci. technol. 4 1989
 !only xy screened and effective coupling along z...
   DO ig = 1, sigma_c_st%ngmt
      qxy     = sqrt(tpiba2*((xq_ibk(1)+g(1, ig))**2 + (xq_ibk(2)+g(2, ig))**2))
      wwp1   = (2*pi*rhoi)/(eps0*meff)*(qxy)
      wwp2   = (2*pi*rhoj)/(eps0*meff)*(qxy)
      alpha   = wwp1 + wwp2*exp(-qxy*z)
      if(qxy.gt.eps8) then
!sterne-static polarizability:
         inveps = 1.0 / (eps0 + 4.0d0*pi*rhoi*e2*((1.0/(meff*qxy*kf**2)) - &
&       (1.0 /(qxy**2.0*kf/2.0))*(real(sqrt(dcmplx((qxy/(2.0*kf))-1.0d0))))))
!need to use inkson's relation here: wwq = alpha*(1-inveps(q,0))^{-1}
!using imaginary frequencies:
        wwq2    = alpha/(1.0 - inveps)
        scrcoul(ig,ig)  = -alpha/(w**2 + wwq2)
      else
        scrcoul(ig,ig)  = 0.0d0
      endif
   ENDDO
 else if(screening.eq.3) then
   DO ig = 1, sigma_c_st%ngmt
!Bilayer plasmon model 
      qxy     = sqrt(tpiba2*((xq_ibk(1)+g(1, ig))**2 + (xq_ibk(2)+g(2, ig))**2))
      wwp1  = ((2*pi*rhoi)/(eps0*meff))*qxy
      wwp2  = ((2*pi*rhoj)/(eps0*meff))*qxy
!Still need to put in wwq modes.
      wwq1   = wwp1/(1-inveps)
      wwq2   = wwp2/(1-inveps)
   if(qxy.gt.eps8) then
        inveps = 1.0 / (eps0 + 4.0d0*pi*rhoi*e2*((1.0/(meff*qxy*kf**2)) - &
&      (1.0 /(qxy**2.0*kf/2.0))*(real(sqrt(dcmplx((qxy/(2.0*kf))-1.0d0))))))
!correct expression for two poles (CAREFUL WITH ROUNDING OF EXPONENTIALS!):
      wwpl   = (wwq1 + wwq2)/2.0 + real(sqrt(dcmplx(((wwq1 - wwq2)/2.0)**2 + wwp1*wwp2*exp(-2*qxy*z))))
      wwmi   = (wwq1 + wwq2)/2.0 - real(sqrt(dcmplx(((wwq1 - wwq2)/2.0)**2 + wwp1*wwp2*exp(-2*qxy*z))))
      A      = (wwpl - wwq2 + wwp2*exp(-2*qxy*z))/(wwpl - wwmi)
      scrcoul(ig,ig) = -A/(w**2 + wwpl) - (1-A)/(w**2 + wwmi)
    else
      scrcoul(ig,ig)  = 0.0d0
    endif
   ENDDO
!else if(screening.eq.4) then
!Multilayer plasmon model gives the expression for stern 2-D dielectric response.
!analytic layered electron gas inverse dielectric function from hawyrlak and co-workers.
!bqg=cosh(qg*d)-((2*pi)/(eps0*x))*((eps0*meff)/(kf*x)-sqrt((eps0*meff/(kf*x))-1))*sinh(qxy*d)
!eps^{-1} = sinh(qd)/\sqrt(b^{2}-1)
!     inveps  = sinh(qg*z)/sqrt(b**2 -1)
!# \inveps(q,0) = sinh(qd)/sqrt(b**2-1)
! fqg  = sinh(x*d)/(sqrt(b(x)**2-1))
!#alpha parameter
! aqg     = ((2*pi*n0)/(eps0*meff)*x)/(1-exp(-2*x*d))
! weff2qg = aqg/(1-fqg)
!#finally the static dielectric fxn is
!invepsqg = 1 - aqg/weff2qg
! else
!     STOP
 endif
END SUBROUTINE mod_diel