!-----------------------------------------------------------------------
  ! 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(ig, xq_ibk, w, inveps, screening)
  USE kinds,      ONLY : DP
  USE constants,  ONLY : e2, fpi, RYTOEV, pi, eps8
  USE gvect,      ONLY : ngm, g
  USE cell_base,     ONLY : tpiba2, tpiba, omega, alat, at

IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: xq_ibk(3)
  COMPLEX(DP) :: w, inveps
  REAL(DP)    :: wwp, eps0, q0, wwq, fac, z, xq(3)
  REAL(DP)    :: wwpi2, wwpj2, qxy, meff, invepsqg
  REAL(DP)    :: rhoi, rhoj , wwqi, wwji
  REAL(DP)    :: wwqi2, wwqj2, wwji2, wwq2
  REAL(DP)    :: A, B, wwmi2, wwpl2 
  REAL(DP)    :: epsq, polq
  REAL(DP)    :: qg, qg2, kf2, kf, alpha
  REAL(DP)    :: expqxy, d, x 
  INTEGER     :: screening, ig
!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    = 16.0/RYTOEV ! plasma frequency in Ry
!     eps0   = 7.4         ! 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   = 1.0d0
     rhoj   = 1.0d0
!z is interlayer distance.
     z = 2.0d0
 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)
!(W-v) = (inveps(w) - delta) v
!Wave vector:
     qg2    = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(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 ) )
     inveps =  - wwp**2.d0/(w**2.d0 + wwq**2.d0)

 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...
     qg2     = (g(1,ig)+xq(1))*2 + (g(2,ig)+xq(2))**2
     qg      = sqrt(tpiba2*qg2)

     wwpi2   = (2*pi)/(eps0*meff)*qg
     wwpj2   = (2*pi)/(eps0*meff)*qg

     qxy     = sqrt(tpiba2*((xq(1)+g(1, ig))**2 + (xq(2)+g(2, ig))**2))
     alpha   = wwpi2 + wwpj2*exp(-qxy*z)

!Sterne polarizability.
!    polq    =  ((eps0*meff)/(kf*x) - sqrt((epsq * meff/(kf*qg))-1))
     epsq    = qg2
     polq    = (-2*kf*meff)/(pi*qg)*((epsq*meff)/(kf*x) - real(sqrt((epsq * meff/(kf*qg))-1)))
!vq is screened coulomb, what is epsq?? need to look at dimensions.
!looks like an energy actually.
!     b       = cosh(qxy*d) - vq*polq*sinh(qxy*z) 
!eps^{-1} = sinh(qd)/\sqrt(b^{2}-1)
     inveps  = sinh(qg*z)/sqrt(b**2 -1)
!need to use inkson's relation here... wwq = alpha*(1-inveps(q,0))^{-1}
     wwq2    = alpha/(1 - inveps)
!Skip through frequencies.
     inveps  = alpha/(w**2 - wwq2)
 else if(screening.eq.3) then
!Bilayer plasmon model 
     qg2     = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2
     qg      = sqrt(tpiba2*qg2)
     wwpi2   = ((2*pi*rhoi)/(eps0*meff))*qg
     wwpj2   = ((2*pi*rhoj)/(eps0*meff))*qg
!Still need to put in wwq modes.
!     wwqi2   =
!     wwqj2   =
     expqxy  = exp(-2*qg*z)
!correct expression for two poles:
     wwpl2   = ((wwqi2+wwqj2)/2) + sqrt(((wwqi2-wwqj2)/2)**2 + wwpi2*wwpj2*expqxy)
     wwmi2   = ((wwqi2+wwqj2)/2) - sqrt(((wwqi2-wwqj2)/2)**2 + wwpi2*wwpj2*expqxy)
!some combination to ensure appropriate simplification:
!     A       = 
!     B       = 1 - A
     inveps  =  A/(w**2 - wwpl2) + B/(w**2- wwmi2)
!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)
!# \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