!-----------------------------------------------------------------------
! 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) :: inveps
REAL(DP) :: w
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 = 12.0_dp/RYTOEV ! plasma frequency in Ry
eps0 = 7.4_dp ! static diel constant of MoS2
q0 = 1.9_dp ! characteristic momentum of MoS2 calculated from Resta paper..
!!!!!!
fac = 1.0_dp/(1.0_dp - 1.0_dp/eps0)
!effective mass annd density of electrons in layer i,j.
meff = 1.0_dp
!or should this be the fourier co-efficient of the density at the particular G vector?
rhoi = 1.0_dp
rhoj = 1.0_dp
!z is interlayer distance.
z = 8.5d0
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.0_dp/(1.0_dp-1.0_dp/eps0)
wwq = wwp * sqrt ( fac * (1.0_dp + (qg/eps0/q0)**2.0_dp ) )
!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.0_dp*pi)/(eps0*meff)*qg
wwpj2 = (2.0_dp*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