!
! Copyright (C) 2001-2014 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 slater (rs, ex, vx)
  !-----------------------------------------------------------------------
  !        Slater exchange with alpha=2/3
  !
  USE kinds, ONLY : DP
#ifdef __LIBXC
  use xc_f90_types_m
  use xc_f90_lib_m
#endif
  implicit none
  real(dp), intent(in) :: rs
  real(dp), intent(out):: ex, vx
#ifdef __LIBXC  
  real(dp):: rho 
  real(dp), parameter :: pi34 = 0.6203504908994d0 ! pi34=(3/4pi)^(1/3)
  integer :: func_id = 1  ! Slater Exchange
  integer :: size = 1
  TYPE(xc_f90_pointer_t) :: xc_func
  TYPE(xc_f90_pointer_t) :: xc_info
  
  rho = (pi34/rs)**3
  call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)    
  call xc_f90_lda_exc_vxc(xc_func, size, rho ,ex, vx)  
  call xc_f90_func_end(xc_func)  
#else
  real(dp), parameter  :: f= -0.687247939924714d0, alpha = 2.0d0/3.0d0
  ! f = -9/8*(3/2pi)^(2/3)
  !
  ex = f * alpha / rs
  vx = 4.d0 / 3.d0 * f * alpha / rs
#endif
  !
  return
end subroutine slater
!
!-----------------------------------------------------------------------
subroutine slater1(rs, ex, vx)
  !-----------------------------------------------------------------------
  !        Slater exchange with alpha=1, corresponding to -1.374/r_s Ry
  !        used to recover old results
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ex, vx
  real(DP), parameter  :: f= -0.687247939924714d0, alpha = 1.0d0
  !
  ex = f * alpha / rs
  vx = 4.d0 / 3.d0 * f * alpha / rs
  !
  return
end subroutine slater1
!
!-----------------------------------------------------------------------
subroutine slater_rxc (rs, ex, vx)
  !-----------------------------------------------------------------------
  !        Slater exchange with alpha=2/3 and Relativistic exchange
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi, c_au
  IMPLICIT none
  real (DP):: rs, ex, vx
  !
  real(DP), PARAMETER :: ZERO=0.D0, ONE=1.D0, PFIVE=.5D0, &
       OPF=1.5D0 !, C014=0.014D0
  real (DP):: trd, ftrd, tftm, a0, alp, z, fz, fzp, vxp, xp, &
       beta, sb, alb, c014
  !
  TRD = ONE/3.d0
  FTRD = 4.d0*TRD
  TFTM = 2**FTRD-2.d0
  A0 = (4.d0/(9.d0*PI))**TRD
  C014= 1.0_DP/a0/c_au
  
  !      X-alpha parameter:
  ALP = 2.d0 * TRD
  
  Z = ZERO
  FZ = ZERO
  FZP = ZERO
  
  VXP = -3.d0*ALP/(2.d0*PI*A0*RS)
  XP = 3.d0*VXP/4.d0
  BETA = C014/RS
  SB = SQRT(1.d0+BETA*BETA)
  ALB = LOG(BETA+SB)
  VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
  XP = XP * (ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
  !  VXF = 2**TRD*VXP
  !  EXF = 2**TRD*XP
  VX = VXP
  EX = XP
END SUBROUTINE slater_rxc

!
!-----------------------------------------------------------------------
  subroutine slaterKZK (rs, ex, vx, vol)
  !-----------------------------------------------------------------------
  !        Slater exchange with alpha=2/3, Kwee, Zhang and Krakauer KE
  !        correction
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ex, vx, dL, vol, ga, pi, a0
  real(DP), parameter  ::  a1 = -2.2037d0, &
              a2 = 0.4710d0, a3 = -0.015d0, ry2h = 0.5d0
  real(DP), parameter  :: f= -0.687247939924714d0, alpha = 2.0d0/3.0d0
  ! f = -9/8*(3/2pi)^(2/3)
  !
  pi = 4.d0 * atan(1.d0)
  a0 = f * alpha * 2.d0

  dL = vol**(1.d0/3.d0)
  ga = 0.5d0 * dL *(3.d0 /pi)**(1.d0/3.d0)
  !
  if ( rs .le. ga) then
   ex = a0 / rs + a1 * rs / dL**2.d0 + a2 * rs**2.d0 / dL**3.d0
   vx = (4.d0 * a0 / rs + 2.d0 * a1 * rs / dL**2.d0 + &
        a2 * rs**2.d0 / dL**3.d0 ) / 3.d0
  else
    ex = a0 / ga + a1 * ga / dL**2.d0 + a2 * ga**2.d0 / dL**3.d0 ! solids
    vx = ex
!   ex = a3 * dL**5.d0 / rs**6.d0     ! molecules
!   vx = 3.d0 * ex  
  endif

  ex = ry2h * ex    ! Ry to Hartree
  vx = ry2h * vx
  !
  return
end subroutine slaterKZK
!
!-----------------------------------------------------------------------
subroutine pz (rs, iflag, ec, vc)
  !-----------------------------------------------------------------------
  !     LDA parameterization from Monte Carlo data
  !     iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981)
  !     iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994)
  !
  USE kinds, ONLY : DP
#ifdef __LIBXC
  use xc_f90_types_m
  use xc_f90_lib_m
#endif
  implicit none
  real(dp), intent(in) :: rs
  real(dp), intent(out):: ec, vc
  integer, intent(in)  :: iflag
#ifdef __LIBXC
  real(dp):: rho 
  real(dp), parameter :: pi34 = 0.6203504908994d0 ! pi34=(3/4pi)^(1/3)
  integer :: func_id = 9   ! Perdew & Zunger
    integer :: size = 1
  TYPE(xc_f90_pointer_t) :: xc_func
  TYPE(xc_f90_pointer_t) :: xc_info

  if (iflag.eq.1)  func_id = 9   ! Perdew & Zunger
  if (iflag.eq.2)  func_id = 11  ! Ortiz & Ballone (PZ)

  rho = (pi34/rs)**3
  call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)    
  call xc_f90_lda_exc_vxc(xc_func, size, rho, ec, vc)  
  call xc_f90_func_end(xc_func)  
#else
  real(DP) :: a (2), b (2), c (2), d (2), gc (2), b1 (2), b2 (2)
  real(DP) :: lnrs, rs12, ox, dox
  !
  data a / 0.0311d0, 0.031091d0 /, b / -0.048d0, -0.046644d0 /, &
       c / 0.0020d0, 0.00419d0 /, d / -0.0116d0, -0.00983d0 /
  data gc / -0.1423d0, -0.103756d0 /, b1 / 1.0529d0, 0.56371d0 /, &
       b2 / 0.3334d0, 0.27358d0 /
  !
  if (rs.lt.1.0d0) then
     ! high density formula
     lnrs = log (rs)
     ec = a (iflag) * lnrs + b (iflag) + c (iflag) * rs * lnrs + d ( &
          iflag) * rs
     vc = a (iflag) * lnrs + (b (iflag) - a (iflag) / 3.d0) + 2.d0 / &
          3.d0 * c (iflag) * rs * lnrs + (2.d0 * d (iflag) - c (iflag) ) &
          / 3.d0 * rs
  else
     ! interpolation formula
     rs12 = sqrt (rs)
     ox = 1.d0 + b1 (iflag) * rs12 + b2 (iflag) * rs
     dox = 1.d0 + 7.d0 / 6.d0 * b1 (iflag) * rs12 + 4.d0 / 3.d0 * &
          b2 (iflag) * rs
     ec = gc (iflag) / ox
     vc = ec * dox / ox
  endif
#endif
  !
  return
end subroutine pz
!
!-----------------------------------------------------------------------
subroutine pzKZK (rs, ec, vc, vol)
  !-----------------------------------------------------------------------
  !     LDA parameterization from Monte Carlo data
  !     iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981)
  !     iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc, ec0 (2), vc0(2), ec0p 
  integer :: iflag, kr
  !
  real(DP) :: a (2), b (2), c (2), d (2), gc (2), b1 (2), b2 (2)
  real(DP) :: lnrs, rs12, ox, dox, lnrsk, rsk
  real(DP) :: a1, grs, g1, g2, g3, g4, dL, vol, gh, gl, grsp
  real(DP) :: f3, f2, f1, f0, pi
  real(DP) :: D1, D2, D3, P1, P2, ry2h
  !
  data a / 0.0311d0, 0.031091d0 /, b / -0.048d0, -0.046644d0 /, &
       c / 0.0020d0, 0.00419d0 /, d / -0.0116d0, -0.00983d0 /
  data gc / -0.1423d0, -0.103756d0 /, b1 / 1.0529d0, 0.56371d0 /, &
       b2 / 0.3334d0, 0.27358d0 /
  data a1 / -2.2037 /, g1 / 0.1182 /, g2 / 1.1656 /, g3 / -5.2884 /, &
       g4 / -1.1233 /
  data ry2h / 0.5d0 /
  !
  iflag = 1
  pi = 4.d0 * atan(1.d0)
  dL = vol**(1.d0/3.d0)
  gh = 0.5d0 * dL / (2.d0 * pi)**(1.d0/3.d0)
  gl = dL * (3.d0 / 2.d0 / pi)**(1.d0/3.d0)

  rsk = gh
  do kr = 1, 2
  lnrsk = log (rsk)
  if (rsk.lt.1.0d0) then
     ! high density formula
    ec0(kr) = a(iflag) *lnrsk + b(iflag) + c(iflag) * rsk * lnrsk + d( &
          iflag) * rsk
    vc0(kr) = a(iflag) * lnrsk + (b(iflag) - a(iflag) / 3.d0) + 2.d0 / &
      3.d0 * c (iflag) * rsk * lnrsk + (2.d0 * d (iflag) - c (iflag) ) &
          / 3.d0 * rsk
  else
     ! interpolation formula
     rs12 = sqrt (rsk)
     ox = 1.d0 + b1 (iflag) * rs12 + b2 (iflag) * rsk
     dox = 1.d0 + 7.d0 / 6.d0 * b1 (iflag) * rs12 + 4.d0 / 3.d0 * &
          b2 (iflag) * rsk
     ec0(kr) = gc (iflag) / ox
     vc0(kr) = ec0(kr) * dox / ox
  endif
  !
   grs  = g1 * rsk * lnrsk + g2 * rsk + g3 * rsk**1.5d0 + g4 * rsk**2.d0
   grsp = g1 * lnrsk + g1 + g2 + 1.5d0 * g3 * rsk**0.5d0 + & 
        2.d0 * g4 * rsk
   ec0(kr)  = ec0(kr) + (-a1 * rsk / dL**2.d0 + grs / dL**3.d0) * ry2h
   vc0(kr)  = vc0(kr) + (-2.d0 * a1 * rsk / dL**2.d0 / 3.d0 + &
           grs / dL**3.d0 -  grsp * rsk / 3.d0 / dL**3.d0) * ry2h
  !
  rsk = rs
  enddo

  lnrs = log (rs)
  if (rs .le. gh) then
   ec = ec0(2)
   vc = vc0(2)
  else
     if ( rs .le. gl) then
        ec0p = 3.d0 * (ec0(1) - vc0(1)) / gh
        P1 = 3.d0 *  ec0(1) - gh * ec0p
        P2 = ec0p
        D1 = gl - gh
        D2 = gl**2.d0 - gh**2.d0
        D3 = gl**3.d0 - gh**3.d0
        f2 = 2.d0 * gl**2.d0 * P2 * D1 + D2 * P1
        f2 = f2 / (-(2.d0*gl*D1)**2.d0 + 4.d0*gl*D1*D2 - D2**2.d0 )
        f3 = - (P2 + 2.d0*D1*f2) / (3.d0 * D2)
        f1 = - (P1 + D2 * f2) / (2.d0 * D1)
        f0 = - gl * (gl * f2 + 2.d0 * f1) / 3.d0
        !
        ec = f3 * rs**3.d0 + f2 * rs**2.d0 + f1 * rs + f0
        vc = f2 * rs**2.d0 / 3.d0 + f1 * 2.d0 * rs / 3.d0 + f0
     else
        ec = 0.d0
        vc = 0.d0
     endif
    endif
  !
  return
end subroutine pzKZK
!
!-----------------------------------------------------------------------
subroutine vwn (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, b, c, x0
  parameter (a = 0.0310907d0, b = 3.72744d0, c = 12.9352d0, x0 = -0.10498d0)
  real(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt
  !
  q = sqrt (4.d0 * c - b * b)
  f1 = 2.d0 * b / q
  f2 = b * x0 / (x0 * x0 + b * x0 + c)
  f3 = 2.d0 * (2.d0 * x0 + b) / q
  rs12 = sqrt (rs)
  fx = rs + b * rs12 + c
  qx = atan (q / (2.d0 * rs12 + b) )
  ec = a * (log (rs / fx) + f1 * qx - f2 * (log ( (rs12 - x0) **2 / &
       fx) + f3 * qx) )
  tx = 2.d0 * rs12 + b
  tt = tx * tx + q * q
  vc = ec - rs12 * a / 6.d0 * (2.d0 / rs12 - tx / fx - 4.d0 * b / &
       tt - f2 * (2.d0 / (rs12 - x0) - tx / fx - 4.d0 * (2.d0 * x0 + b) &
       / tt) )
  !
  return
end subroutine vwn

!-----------------------------------------------------------------------
subroutine vwn1_rpa (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, b, c, x0
  parameter (a = 0.0310907_dp, b = 13.0720_dp, c = 42.7198_dp, x0 = -0.409286_dp)
  real(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt
  !
  q = sqrt (4.d0 * c - b * b)
  f1 = 2.d0 * b / q
  f2 = b * x0 / (x0 * x0 + b * x0 + c)
  f3 = 2.d0 * (2.d0 * x0 + b) / q
  rs12 = sqrt (rs)
  fx = rs + b * rs12 + c
  qx = atan (q / (2.d0 * rs12 + b) )
  ec = a * (log (rs / fx) + f1 * qx - f2 * (log ( (rs12 - x0) **2 / &
       fx) + f3 * qx) )
  tx = 2.d0 * rs12 + b
  tt = tx * tx + q * q
  vc = ec - rs12 * a / 6.d0 * (2.d0 / rs12 - tx / fx - 4.d0 * b / &
       tt - f2 * (2.d0 / (rs12 - x0) - tx / fx - 4.d0 * (2.d0 * x0 + b) &
       / tt) )
  !
  return
end subroutine vwn1_rpa

!-----------------------------------------------------------------------
subroutine lyp (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     C. Lee, W. Yang, and R.G. Parr, PRB 37, 785 (1988)
  !     LDA part only
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, b, c, d, pi43
  parameter (a = 0.04918d0, b = 0.132d0 * 2.87123400018819108d0)
  ! pi43 = (4pi/3)^(1/3)
  parameter (pi43 = 1.61199195401647d0, c = 0.2533d0 * pi43, d = &
       0.349d0 * pi43)
  real(DP) :: ecrs, ox
  !
  ecrs = b * exp ( - c * rs)
  ox = 1.d0 / (1.d0 + d * rs)
  ec = - a * ox * (1.d0 + ecrs)
  vc = ec - rs / 3.d0 * a * ox * (d * ox + ecrs * (d * ox + c) )
  !
  return
end subroutine lyp
!
!-----------------------------------------------------------------------
subroutine pw (rs, iflag, ec, vc)
  !-----------------------------------------------------------------------
  !     iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !     iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994)
  !
  USE kinds, ONLY : DP
#ifdef __LIBXC
  use xc_f90_types_m
  use xc_f90_lib_m
#endif
  implicit none
  real(dp), intent(in) :: rs
  real(dp), intent(out):: ec, vc
  integer, intent(in) :: iflag 
#ifdef __LIBXC
  real(dp):: rho 
  real(dp), parameter :: pi34 = 0.6203504908994d0 ! pi34=(3/4pi)^(1/3)
  integer :: func_id = 12   ! Perdew & Zunger
  integer :: size = 1
  TYPE(xc_f90_pointer_t) :: xc_func
  TYPE(xc_f90_pointer_t) :: xc_info

  if (iflag.eq.1)  func_id = 12  ! Perdew & Wang
  if (iflag.eq.2)  func_id = 14  ! Ortiz & Ballone (PW)

  rho =  (pi34/rs)**3
  call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)    
  call xc_f90_lda_exc_vxc(xc_func,size , rho, ec, vc)
  call xc_f90_func_end(xc_func)  
#else
  real(DP) :: a, b1, b2, c0, c1, c2, c3, d0, d1
  parameter (a = 0.031091d0, b1 = 7.5957d0, b2 = 3.5876d0, c0 = a, &
       c1 = 0.046644d0, c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, &
       d1 = 1.4408d0)
  real(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog
  real(DP) :: a1 (2), b3 (2), b4 (2)
  data a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, &
       b4 / 0.49294d0, 0.13354d0 /
  !
  ! high- and low-density formulae implemented but not used in PW case
  ! (reason: inconsistencies in PBE/PW91 functionals)
  !
  if (rs.lt.1d0.and.iflag.eq.2) then
     ! high density formula
     lnrs = log (rs)
     ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs
     vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * &
          lnrs - (2.d0 * c3 + c2) / 3.d0 * rs
  elseif (rs.gt.100.d0.and.iflag.eq.2) then
     ! low density formula
     ec = - d0 / rs + d1 / rs**1.5d0
     vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0
  else
     ! interpolation formula
     rs12 = sqrt (rs)
     rs32 = rs * rs12
     rs2 = rs**2
     om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 (iflag) * rs32 + b4 ( &
          iflag) * rs2)
     dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 ( &
          iflag) * rs32 + 2.d0 * b4 (iflag) * rs2)
     olog = log (1.d0 + 1.0d0 / om)
     ec = - 2.d0 * a * (1.d0 + a1 (iflag) * rs) * olog
     vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 (iflag) * rs) &
          * olog - 2.d0 / 3.d0 * a * (1.d0 + a1 (iflag) * rs) * dom / &
          (om * (om + 1.d0) )
  endif
#endif
  return
end subroutine pw
!
!-----------------------------------------------------------------------
subroutine wigner (rs, ec, vc)
  !-----------------------------------------------------------------------
  !        Wigner correlation
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: pi34, rho13
  parameter (pi34 = 0.6203504908994d0)
  ! pi34=(3/4pi)^(1/3), rho13=rho^(1/3)
  !
  rho13 = pi34 / rs
  vc = - rho13 * ( (0.943656d0 + 8.8963d0 * rho13) / (1.d0 + &
       12.57d0 * rho13) **2)
  ec = - 0.738d0 * rho13 * (0.959d0 / (1.d0 + 12.57d0 * rho13) )
  !
  return
end subroutine wigner
!
!-----------------------------------------------------------------------
subroutine hl (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     L. Hedin and  B.I. Lundqvist,  J. Phys. C 4, 2064 (1971)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, x
  !
  a = log (1.0d0 + 21.d0 / rs)
  x = rs / 21.0d0
  ec = a + (x**3 * a - x * x) + x / 2.d0 - 1.0d0 / 3.0d0
  ec = - 0.0225d0 * ec
  vc = - 0.0225d0 * a
  !
  return
end subroutine hl
!
!-----------------------------------------------------------------------
subroutine gl (rs, ec, vc)
  !-----------------------------------------------------------------------
  !  O. Gunnarsson and B. I. Lundqvist, PRB 13, 4274 (1976)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, vc, ec
  real(DP) :: c, r, x
  parameter (c = 0.0333d0, r = 11.4d0)
  ! c=0.0203, r=15.9 for the paramagnetic case
  !
  x = rs / r
  vc = - c * log (1.d0 + 1.d0 / x)
  ec = - c * ( (1.d0 + x**3) * log (1.d0 + 1.d0 / x) - 1.0d0 / &
       3.0d0 + x * (0.5d0 - x) )
  !
  return
end subroutine gl
!
!-----------------------------------------------------------------------
subroutine becke86b(rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! Becke 1986 gradient correction to exchange
  ! A.D. Becke, J. Chem. Phys. 85 (1986) 7184
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  real(DP) :: arho, agrho, beta, gamma
  parameter (beta = 0.00375_dp, gamma=0.007_dp)
  real(dp) :: sgp1, sgp1_45, sgp1_95
  real(dp) :: rdg2_43, rdg2_73, rdg2_83, rdg2_4, rdg4_5

  arho = 0.5d0 * rho 
  agrho = 0.25d0 * grho 

  rdg2_43 = agrho / arho**(4d0/3d0)
  rdg2_73 = rdg2_43 / arho
  rdg2_83 = rdg2_43 * rdg2_43 / agrho
  rdg2_4 = rdg2_43 * rdg2_83 / agrho
  rdg4_5 = rdg2_73 * rdg2_83

  sgp1 = 1d0 + gamma * rdg2_83
  sgp1_45 = sgp1**(-4d0/5d0)
  sgp1_95 = sgp1_45 / sgp1

  sx = -2d0 * beta * agrho / arho**(4d0/3d0) * sgp1_45
  v1x = -beta * (-4d0/3d0 * rdg2_73 * sgp1_45 + 32d0/15d0 * gamma * rdg4_5 * sgp1_95)
  v2x = -beta * (sgp1_45 * rdg2_43 / agrho - 4d0/5d0 * gamma * rdg2_4 * sgp1_95)

end subroutine becke86b
!
!-----------------------------------------------------------------------
subroutine becke88 (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! Becke exchange: A.D. Becke, PRA 38, 3098 (1988)
  ! only gradient-corrected part, no Slater term included
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  real(DP) :: beta, third, two13
  parameter (beta = 0.0042d0)
  parameter (third = 1.d0 / 3.d0, two13 = 1.259921049894873d0)
  ! two13 = 2^(1/3)
  real(DP) :: rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee
  !
  rho13 = rho**third
  rho43 = rho13**4
  xs = two13 * sqrt (grho) / rho43
  xs2 = xs * xs
  sa2b8 = sqrt (1.0d0 + xs2)
  shm1 = log (xs + sa2b8)
  dd = 1.0d0 + 6.0d0 * beta * xs * shm1
  dd2 = dd * dd
  ee = 6.0d0 * beta * xs2 / sa2b8 - 1.d0
  sx = two13 * grho / rho43 * ( - beta / dd)
  v1x = - (4.d0 / 3.d0) / two13 * xs2 * beta * rho13 * ee / dd2
  v2x = two13 * beta * (ee-dd) / (rho43 * dd2)
  !
  return
end subroutine becke88
!
!-----------------------------------------------------------------------
subroutine ggax (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! Perdew-Wang GGA (PW91), exchange part:
  ! J.P. Perdew et al.,PRB 46, 6671 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  real(DP) :: f1, f2, f3, f4, f5
  parameter (f1 = 0.19645d0, f2 = 7.7956d0, f3 = 0.2743d0, f4 = &
       0.1508d0, f5 = 0.004d0)
  real(DP) :: fp1, fp2
  parameter (fp1 = -0.019292021296426d0, fp2 = 0.161620459673995d0)
  ! fp1 = -3/(16 pi)*(3 pi^2)^(-1/3)
  ! fp2 = (1/2)(3 pi^2)**(-1/3)
  real(DP) :: rhom43, s, s2, s3, s4, exps, as, sa2b8, shm1, bs, das, &
       dbs, dls
  !
  rhom43 = rho** ( - 4.d0 / 3.d0)
  s = fp2 * sqrt (grho) * rhom43
  s2 = s * s
  s3 = s2 * s
  s4 = s2 * s2
  exps = f4 * exp ( - 100.d0 * s2)
  as = f3 - exps - f5 * s2
  sa2b8 = sqrt (1.0d0 + f2 * f2 * s2)
  shm1 = log (f2 * s + sa2b8)
  bs = 1.d0 + f1 * s * shm1 + f5 * s4
  das = (200.d0 * exps - 2.d0 * f5) * s
  dbs = f1 * (shm1 + f2 * s / sa2b8) + 4.d0 * f5 * s3
  dls = (das / as - dbs / bs)
  sx = fp1 * grho * rhom43 * as / bs
  v1x = - 4.d0 / 3.d0 * sx / rho * (1.d0 + s * dls)
  v2x = fp1 * rhom43 * as / bs * (2.d0 + s * dls)
  !
  return
end subroutine ggax

!
!-----------------------------------------------------------------------
subroutine rPW86 (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! PRB 33, 8800 (1986) and J. Chem. Theory comp. 5, 2754 (2009)
  !
  USE kinds
  implicit none

  real(DP), intent(in) :: rho, grho
  real(DP), intent(out) :: sx, v1x, v2x
  real(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, grad_rho, df_ds
  real(DP) :: a, b, c, s_prefactor, Ax, four_thirds
  parameter( a = 1.851d0, b = 17.33d0, c = 0.163d0, s_prefactor = 6.18733545256027d0, &
       Ax = -0.738558766382022d0, four_thirds = 4.d0/3.d0)

  grad_rho = sqrt(grho)

  s = grad_rho/(s_prefactor*rho**(four_thirds))

  s_2 = s**2
  s_3 = s_2 * s
  s_4 = s_2**2
  s_5 = s_3 * s_2
  s_6 = s_2 * s_4

  !! Calculation of energy
  fs = (1 + a*s_2 + b*s_4 + c*s_6)**(1.d0/15.d0)
  sx = Ax * rho**(four_thirds) * (fs -1.0D0)


  !! Calculation of the potential
  df_ds = (1.d0/(15.d0*fs**(14.0D0)))*(2*a*s + 4*b*s_3 + 6*c*s_5)


  v1x = Ax*(four_thirds)*(rho**(1.d0/3.d0)*(fs -1.0D0) &
       -grad_rho/(s_prefactor * rho)*df_ds)

  v2x = Ax * df_ds/(s_prefactor*grad_rho)

end subroutine rPW86

subroutine PW86 (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! Perdew-Wang 1986 exchange gradient correction: PRB 33, 8800 (1986) 
  ! 
  USE kinds
  implicit none

  real(DP), intent(in) :: rho, grho
  real(DP), intent(out) :: sx, v1x, v2x
  real(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, grad_rho, df_ds
  real(DP) :: a, b, c, s_prefactor, Ax, four_thirds
  parameter( a = 1.296_dp, b = 14_dp, c = 0.2_dp, s_prefactor = 6.18733545256027_dp, &
       Ax = -0.738558766382022_dp, four_thirds = 4._dp/3._dp)

  grad_rho = sqrt(grho)

  s = grad_rho/(s_prefactor*rho**(four_thirds))

  s_2 = s**2
  s_3 = s_2 * s
  s_4 = s_2**2
  s_5 = s_3 * s_2
  s_6 = s_2 * s_4

  !! Calculation of energy
  fs = (1 + a*s_2 + b*s_4 + c*s_6)**(1.d0/15.d0)
  sx = Ax * rho**(four_thirds) * (fs-1d0)

  !! Calculation of the potential
  df_ds = (1.d0/(15.d0*fs**(14d0)))*(2*a*s + 4*b*s_3 + 6*c*s_5)

  v1x = Ax*(four_thirds)*(rho**(1.d0/3.d0)*(fs-1d0) &
       -grad_rho/(s_prefactor * rho)*df_ds)

  v2x = Ax * df_ds/(s_prefactor*grad_rho)

end subroutine PW86
!
!-----------------------------------------------------------------------
subroutine cx13 (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! The new exchange partner for a vdW-DF1-cx suggested
  ! by K. Berland and P. Hyldgaard, see PRB 89, 035412 (2014), 
  ! to test the plasmon nature of the vdW-DF1 inner functional.
  !
  USE kinds
  implicit none
  real(DP), intent(in) :: rho, grho
  real(DP), intent(out) :: sx, v1x, v2x
  real(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, fs_rPW86, df_rPW86_ds, grad_rho, df_ds
  real(DP) :: alp, beta, a, b, c, s_prefactor, Ax, four_thirds, mu_LM
  parameter( alp = 0.021789d0, beta=1.15d0,  a = 1.851d0, b = 17.33d0, c = 0.163d0, mu_LM = 0.09434d0, &
       s_prefactor = 6.18733545256027d0, &
       Ax = -0.738558766382022d0, four_thirds = 4.d0/3.d0)

  grad_rho = sqrt(grho)

  s = grad_rho/(s_prefactor*rho**(four_thirds))

  s_2 = s*s
  s_3 = s_2 * s
  s_4 = s_2 * s_2
  s_5 = s_3 * s_2
  s_6 = s_2 * s_2 *s_2


  !! Calculation of energy
  fs_rPW86 = (1 + a*s_2 + b*s_4 + c*s_6)**(1.d0/15.d0)
  fs = 1.d0/(1+alp*s_6) * (1 + mu_LM *s_2) &
       +  alp*s_6/(beta+alp*s_6)*fs_rPW86

  sx = Ax * rho**(four_thirds) * (fs -1.0D0)


  !! Calculation of the potential
  df_rPW86_ds = (1.d0/(15.d0*fs_rPW86**(14.0D0)))*(2*a*s + 4*b*s_3 + 6*c*s_5)

  df_ds = 1.d0/(1+alp*s_6)**2*( 2.0d0*mu_LM*s*(1+alp*s_6) - 6.0d0*alp*s_5*( 1 + mu_LM*s_2))  &
         + alp*s_6/(beta+alp*s_6)*df_rPW86_ds  &
         + 6.0d0*alp*s_5*fs_rPW86/(beta+alp*s_6)*(1- alp*s_6/(beta + alp*s_6))

  v1x = Ax*(four_thirds)*(rho**(1.d0/3.d0)*(fs -1.0D0) &
       -grad_rho/(s_prefactor * rho)*df_ds)
  v2x = Ax * df_ds/(s_prefactor*grad_rho)

end subroutine cx13
!
!---------------------------------------------------------------
subroutine c09x (rho, grho, sx, v1x, v2x)
  !---------------------------------------------------------------
  ! Cooper '09 exchange for vdW-DF (without Slater exchange):
  ! V. R. Cooper, Phys. Rev. B 81, 161104(R) (2010)
  !
  ! Developed thanks to the contribution of 
  ! Ikutaro Hamada - ikutaro@wpi-aimr.tohoku.ac.jp
  ! WPI-Advanced Institute of Materials Research, Tohoku University  
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  ! input: charge and squared gradient
  ! output: energy
  ! output: potential
  ! local variables
  real(DP) :: kf, agrho, s1, s2, ds, dsg, exunif, fx
  ! (3*pi2*|rho|)^(1/3)
  ! |grho|
  ! |grho|/(2*kf*|rho|)
  ! s^2
  ! n*ds/dn
  ! n*ds/d(gn)
  ! exchange energy LDA part
  ! exchange energy gradient part
  real(DP) :: dxunif, dfx, f1, f2, f3, dfx1, dfx2
  ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
  real(DP) :: third, c1, c2, c5
  parameter (third = 1.d0 / 3.d0, c1 = 0.75d0 / pi , &
       c2 = 3.093667726280136d0, c5 = 4.d0 * third)
  ! parameters of the functional
  real(DP) :: kappa, mu, alpha
  data kappa / 1.245d0  /, &
       mu  / 0.0617d0 /, &
       alpha / 0.0483d0 /
  !
  agrho = sqrt (grho)
  kf = c2 * rho**third
  dsg = 0.5d0 / kf
  s1 = agrho * dsg / rho
  s2 = s1 * s1
  ds = - c5 * s1
  !
  !   Energy
  !
  f1 = exp( - alpha * s2 )
  f2 = exp( - alpha * s2 / 2.0d0 )
  f3 = mu * s2 * f1
  fx = f3 + kappa * ( 1.0d0 - f2 )
  exunif = - c1 * kf
  sx = exunif * fx
  !
  !   Potential
  !
  dxunif = exunif * third
  dfx1 = 2.0d0 * mu * s1 * ( 1.0d0 - alpha * s2 ) * f1
  dfx2 = kappa * alpha * s1 * f2 
  dfx = dfx1 + dfx2 
  v1x = sx + dxunif * fx + exunif * dfx * ds
  v2x = exunif * dfx * dsg / agrho

  sx = sx * rho
  return
end subroutine c09x
!---------------------------------------------------------------
subroutine b86b (rho, grho, iflag, sx, v1x, v2x)
  !---------------------------------------------------------------
  ! Becke exchange (without Slater exchange):
  ! iflag=1: A. D. Becke, J. Chem. Phys. 85, 7184 (1986) (B86b)
  ! iflag=2: J. Klimes, Phys. Rev. B 83, 195131 (2011). (OptB86b)
  ! iflag=3: I. Hamada, Phys. Rev. B 89, 121103(R) (B86R)
  !
  ! Ikutaro Hamada - HAMADA.Ikutaro@nims.go.jp
  ! National Institute for Materials Science
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  ! input: charge and squared gradient
  ! output: energy
  ! output: potential
  integer :: iflag
  ! local variables
  real(DP) :: kf, agrho, s1, s2, ds, dsg, exunif, fx
  ! (3*pi2*|rho|)^(1/3)
  ! |grho|
  ! |grho|/(2*kf*|rho|)
  ! s^2
  ! n*ds/dn
  ! n*ds/d(gn)
  ! exchange energy LDA part
  ! exchange energy gradient part
  real(DP) :: dxunif, dfx, f1, f2, f3, dfx1
  ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
  real(DP) :: third, c1, c2, c5
  parameter (third = 1.d0 / 3.d0, c1 = 0.75d0 / pi , &
       c2 = 3.093667726280136d0, c5 = 4.d0 * third)
  ! parameters of the functional
  real(DP) :: k (3), mu(3)
  data k / 0.5757d0, 1.0000d0, 0.711357d0/, &
       mu/ 0.2449d0, 0.1234d0, 0.1234d0  /
  !
  agrho = sqrt (grho)
  kf = c2 * rho**third
  dsg = 0.5d0 / kf
  s1 = agrho * dsg / rho
  s2 = s1 * s1
  ds = - c5 * s1
  !
  !   Energy
  !
  f1=mu(iflag)*s2
  f2=1.d0+mu(iflag)*s2/k(iflag)
  f3=f2**(4.d0/5.d0)
  fx= f1/f3
  exunif = - c1 * kf
  sx = exunif * fx
  !
  !   Potential
  !
  dxunif = exunif * third
  dfx1 = 1.d0+(1.d0/5.d0)*mu(iflag)*s2/k(iflag)
  dfx = 2.d0* mu(iflag) * s1 * dfx1 / (f2 * f3)
  v1x = sx + dxunif * fx + exunif * dfx * ds
  v2x = exunif * dfx * dsg / agrho
  sx = sx * rho
  return
end subroutine b86b
!
!-----------------------------------------------------------------------
subroutine perdew86 (rho, grho, sc, v1c, v2c)
  !-----------------------------------------------------------------------
  ! Perdew gradient correction on correlation: PRB 33, 8822 (1986)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sc, v1c, v2c
  real(DP) :: p1, p2, p3, p4, pc1, pc2, pci
  parameter (p1 = 0.023266d0, p2 = 7.389d-6, p3 = 8.723d0, p4 = &
       0.472d0)
  parameter (pc1 = 0.001667d0, pc2 = 0.002568d0, pci = pc1 + pc2)
  real(DP) :: third, pi34
  parameter (third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0)
  ! pi34=(3/4pi)^(1/3)
  real(DP) :: rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs
  real(DP) :: dcna, dcnb, dcn, phi, ephi
  !
  rho13 = rho**third
  rho43 = rho13**4
  rs = pi34 / rho13
  rs2 = rs * rs
  rs3 = rs * rs2
  cna = pc2 + p1 * rs + p2 * rs2
  cnb = 1.d0 + p3 * rs + p4 * rs2 + 1.d4 * p2 * rs3
  cn = pc1 + cna / cnb
  drs = - third * pi34 / rho43
  dcna = (p1 + 2.d0 * p2 * rs) * drs
  dcnb = (p3 + 2.d0 * p4 * rs + 3.d4 * p2 * rs2) * drs
  dcn = dcna / cnb - cna / (cnb * cnb) * dcnb
  phi = 0.192d0 * pci / cn * sqrt (grho) * rho** ( - 7.d0 / 6.d0)
  ! SdG: in the original paper 1.745*0.11=0.19195 is used
  ephi = exp ( - phi)
  sc = grho / rho43 * cn * ephi
  v1c = sc * ( (1.d0 + phi) * dcn / cn - ( (4.d0 / 3.d0) - (7.d0 / &
       6.d0) * phi) / rho)
  v2c = cn * ephi / rho43 * (2.d0 - phi)
  !
  return
end subroutine perdew86
!
!-----------------------------------------------------------------------
subroutine glyp (rho, grho, sc, v1c, v2c)
  !-----------------------------------------------------------------------
  ! Lee Yang Parr: gradient correction part
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sc, v1c, v2c
  real(DP) :: a, b, c, d
  parameter (a = 0.04918d0, b = 0.132d0, c = 0.2533d0, d = 0.349d0)
  real(DP) :: rhom13, rhom43, rhom53, om, xl, ff, dom, dxl
  !
  rhom13 = rho** ( - 1.d0 / 3.d0)
  om = exp ( - c * rhom13) / (1.d0 + d * rhom13)
  xl = 1.d0 + (7.d0 / 3.d0) * (c * rhom13 + d * rhom13 / (1.d0 + d * &
       rhom13) )
  ff = a * b * grho / 24.d0
  rhom53 = rhom13**5
  sc = ff * rhom53 * om * xl
  dom = - om * (c + d+c * d * rhom13) / (1.d0 + d * rhom13)
  dxl = (7.d0 / 3.d0) * (c + d+2.d0 * c * d * rhom13 + c * d * d * &
       rhom13**2) / (1.d0 + d * rhom13) **2
  rhom43 = rhom13**4
  v1c = - ff * rhom43 / 3.d0 * (5.d0 * rhom43 * om * xl + rhom53 * &
       dom * xl + rhom53 * om * dxl)
  v2c = 2.d0 * sc / grho
  !
  return
end subroutine glyp
!
!-----------------------------------------------------------------------
subroutine ggac (rho, grho, sc, v1c, v2c)
  !-----------------------------------------------------------------------
  ! Perdew-Wang GGA (PW91) correlation part
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sc, v1c, v2c
  real(DP) :: al, pa, pb, pc, pd, cx, cxc0, cc0
  parameter (al = 0.09d0, pa = 0.023266d0, pb = 7.389d-6, pc = &
       8.723d0, pd = 0.472d0)
  parameter (cx = -0.001667d0, cxc0 = 0.002568d0, cc0 = - cx + cxc0)
  real(DP) :: third, pi34, nu, be, xkf, xks
  parameter (third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0)
  parameter (nu = 15.755920349483144d0, be = nu * cc0)
  parameter (xkf = 1.919158292677513d0, xks = 1.128379167095513d0)
  ! pi34=(3/4pi)^(1/3),  nu=(16/pi)*(3 pi^2)^(1/3)
  ! xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
  real(DP) :: kf, ks, rs, rs2, rs3, ec, vc, t, expe, af, bf, y, xy, &
       qy, s1
  real(DP) :: h0, dh0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, &
       dh1, ddh1
  !
  rs = pi34 / rho**third
  rs2 = rs * rs
  rs3 = rs * rs2
  call pw (rs, 1, ec, vc)
  kf = xkf / rs
  ks = xks * sqrt (kf)
  t = sqrt (grho) / (2.d0 * ks * rho)
  expe = exp ( - 2.d0 * al * ec / (be * be) )
  af = 2.d0 * al / be * (1.d0 / (expe-1.d0) )
  bf = expe * (vc - ec)
  y = af * t * t
  xy = (1.d0 + y) / (1.d0 + y + y * y)
  qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
  s1 = 1.d0 + 2.d0 * al / be * t * t * xy
  h0 = be * be / (2.d0 * al) * log (s1)
  dh0 = be * t * t / s1 * ( - 7.d0 / 3.d0 * xy - qy * (af * bf / &
       be-7.d0 / 3.d0) )
  ddh0 = be / (2.d0 * ks * ks * rho) * (xy - qy) / s1
  ee = - 100.d0 * (ks / kf * t) **2
  cna = cxc0 + pa * rs + pb * rs2
  dcna = pa * rs + 2.d0 * pb * rs2
  cnb = 1.d0 + pc * rs + pd * rs2 + 1.d4 * pb * rs3
  dcnb = pc * rs + 2.d0 * pd * rs2 + 3.d4 * pb * rs3
  cn = cna / cnb - cx
  dcn = dcna / cnb - cna * dcnb / (cnb * cnb)
  h1 = nu * (cn - cc0 - 3.d0 / 7.d0 * cx) * t * t * exp (ee)
  dh1 = - third * (h1 * (7.d0 + 8.d0 * ee) + nu * t * t * exp (ee) &
       * dcn)
  ddh1 = 2.d0 * h1 * (1.d0 + ee) * rho / grho
  sc = rho * (h0 + h1)
  v1c = h0 + h1 + dh0 + dh1
  v2c = ddh0 + ddh1
  !
  return
end subroutine ggac
!
!---------------------------------------------------------------
subroutine pbex (rho, grho, iflag, sx, v1x, v2x)
  !---------------------------------------------------------------
  !
  ! PBE exchange (without Slater exchange):
  ! iflag=1  J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)
  ! iflag=2  "revised' PBE: Y. Zhang et al., PRL 80, 890 (1998)
  ! iflag=3  PBEsol: J.P.Perdew et al., PRL 100, 136406 (2008)
  ! iflag=4  PBEQ2D: L. Chiodo et al., PRL 108, 126402 (2012)
  ! iflag=5  optB88: Klimes et al., J. Phys. Cond. Matter, 22, 022201 (2010)
  ! iflag=6  optB86b: Klimes et al., Phys. Rev. B 83, 195131 (2011)
  ! iflag=7  ev: Engel and Vosko, PRB 47, 13164 (1991)
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
#ifdef __LIBXC
  use xc_f90_types_m
  use xc_f90_lib_m
#endif
  implicit none
  real(dp), intent(in) :: rho, grho
  ! input: charge and squared gradient
  real(dp), intent(out):: sx, v1x, v2x
  ! output: energy, potential
  integer, intent(in) :: iflag
#ifdef __LIBXC
  ! local variables
  integer :: func_id = -1 ! not set
  integer :: size = 1
  TYPE(xc_f90_pointer_t) :: xc_func
  TYPE(xc_f90_pointer_t) :: xc_info
  real(dp) :: exc, ex_lda = 0.0d0 , vx_lda = 0.0d0

  if (iflag.eq.1)  func_id = 101
  if (iflag.eq.2)  func_id = 102
  if (iflag.eq.3)  func_id = 116
  if (iflag.eq.5)  func_id = 141
  if (func_id==-1) call errore('pbex','case not implemented with libxc',iflag)

  call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)    
  call xc_f90_gga_exc_vxc(xc_func, size, rho, grho, exc, v1x, v2x)  
  call xc_f90_func_end(xc_func)  

  ! remove Slater term for compatibility with QE  
  call xc_f90_func_init(xc_func, xc_info, 1, XC_UNPOLARIZED)       
  call xc_f90_lda_exc_vxc(xc_func, size, rho, ex_lda, vx_lda)  
  call xc_f90_func_end(xc_func) 
  exc = exc - ex_lda 
  v1x = v1x - vx_lda 
  
  sx = exc * rho  ! e_x = rho * \epsilon_x
  v2x = v2x*2.0_dp

#else
  ! local variables
  real(DP) :: kf, agrho, s1, s2, ds, dsg, exunif, fx
  ! (3*pi2*|rho|)^(1/3)
  ! |grho|
  ! |grho|/(2*kf*|rho|)
  ! s^2
  ! n*ds/dn
  ! n*ds/d(gn)
  ! exchange energy LDA part
  ! exchange energy gradient part
  real(DP) :: dxunif, dfx, f1, f2, f3, dfx1
  real(DP) :: p,  amu, ab, c, dfxdp, dfxds, upbe, uge, s, ak, aa
  ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
  real(DP), parameter :: third = 1._DP / 3._DP, c1 = 0.75_DP / pi , &
       c2 = 3.093667726280136_DP, c5 = 4._DP * third, &
       c6 = c2*2.51984210, c7=5._DP/6._DP, c8=0.8_DP ! (3pi^2)^(1/3)*2^(4/3)
  ! parameters of the functional
  real(DP) :: k (6), mu(6), ev(6)
  !           pbe        rpbe        pbesol   pbeq2d      optB88  optB86b
  data k / 0.804_DP,   1.2450D0,   0.804_DP , 0.804_DP ,    0.0 ,  0.0 /, &
       mu/ 0.21951_DP, 0.21951_DP, 0.12345679012345679012_DP,             &
                                   0.12345679012345679,     0.22 , 0.1234/, &
       ev / 1.647127_DP, 0.980118_DP, 0.017399_DP, 1.523671_DP, 0.367229_DP, &
                                   0.011282_DP /  ! a and b parameters of Engel and Vosko
  agrho = sqrt (grho)
  kf = c2 * rho**third
  dsg = 0.5_DP / kf
  s1 = agrho * dsg / rho
  s2 = s1 * s1
  ds = - c5 * s1
  !
  !   Energy
  !
  if ( iflag == 4) then
     p=s1*s1
     s=s1
     ak=0.804_DP
     amu=10._DP/81._DP
     ab=0.5217_DP
     c=2._DP
     fx =  ak - ak / (1.0_dp + amu * p / ak)  + p**2 * (1 + p) &
           /(10**c + p**3) * (-1.0_dp - ak + ak / (1.0_dp + amu * p / ak) &
           + ab * p ** (-0.1d1/ 0.4D1))
  elseif ( iflag == 5) then
     ab=mu(iflag)*c7 ! mu/ab=1.2
     p=s1*c6
     c=log(p+sqrt(p*p+1)) ! asinh(p)
     dfx1=1+ab*s1*c
     fx =  mu(iflag)*s1*s1/dfx1
  elseif ( iflag == 6) then
     p=mu(iflag)*s1*s1
     fx =  p / ( 1 + p )**c8
  elseif ( iflag == 7) then
     s=s2*s2
     f1 =  1 + ev(1)*s2 + ev(2)*s + ev(3)*s*s2
     f2 =  1 + ev(4)*s2 + ev(5)*s + ev(6)*s*s2
     fx = f1 / f2 - 1
  else
     f1 = s2 * mu(iflag) / k (iflag)
     f2 = 1._DP + f1
     f3 = k (iflag) / f2
     fx = k (iflag) - f3
  end if
  exunif = - c1 * kf
  sx = exunif * fx
  !
  !   Potential
  !
  dxunif = exunif * third
  if ( iflag == 4) then
      dfxdp = dble(1 / (1 + amu * p / ak) ** 2 * amu) + dble(2 * p * (1 &
     + p) / (10 ** c + p ** 3) * (-1 - ak + ak / (1 + amu * p / ak) + ab &
      * p ** (-0.1d1 / 0.4D1))) + dble(p ** 2 / (10 ** c + p ** 3) * ( &
     -1 - ak + ak / (1 + amu * p / ak) + ab * p ** (-0.1d1 / 0.4D1))) - &
      dble(3 * p ** 4 * (1 + p) / (10 ** c + p ** 3) ** 2 * (-1 - ak + &
     ak / (1 + amu * p / ak) + ab * p ** (-0.1d1 / 0.4D1))) + dble(p ** &
      2) * dble(1 + p) / dble(10 ** c + p ** 3) * (-dble(1 / (1 + amu * &
      p / ak) ** 2 * amu) - dble(ab * p ** (-0.5d1 / 0.4D1)) / 0.4D1)

      dfxds=dfxdp*2._DP*s
      dfx=dfxds
  elseif (iflag == 5) then
     dfx=2*fx/s1-fx/dfx1*(ab*c+ab*s1/sqrt(p*p+1)*c6)
  elseif (iflag == 6) then
     dfx=2*mu(iflag)*s1*fx*(1+(1-c8)*p)/(p*(1+p))
  elseif (iflag == 7) then
    dfx  =  ev(1) + 2*ev(2)*s2 + 3*ev(3)*s  
    dfx1 =  ev(4) + 2*ev(5)*s2 + 3*ev(6)*s 
    dfx  = 2 * s1 * ( dfx - f1*dfx1/f2 ) / f2
  else
     dfx1 = f2 * f2
     dfx = 2._DP * mu(iflag) * s1 / dfx1
  end if
  v1x = sx + dxunif * fx + exunif * dfx * ds
  v2x = exunif * dfx * dsg / agrho
  sx = sx * rho
#endif
  return
end subroutine pbex
!
!---------------------------------------------------------------
subroutine pbex_vec (rho, grho, iflag, sx, v1x, v2x, length, small)
  !---------------------------------------------------------------
  !
  ! PBE exchange (without Slater exchange):
  ! iflag=1  J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)
  ! iflag=2  "revised' PBE: Y. Zhang et al., PRL 80, 890 (1998)
  ! iflag=3  PBEsol: J.P.Perdew et al., PRL 100, 136406 (2008)
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
  implicit none
  integer, intent(in)   :: length
  integer, intent(in)   :: iflag
  real(DP), intent(in)  :: small
  real(DP), intent(in)  :: rho(length), grho(length)
  real(DP), intent(out) :: sx(length), v1x(length), v2x(length)
  ! local variables
  integer :: i
  real(DP) :: kf, agrho, s1, dsg, exunif, fx
  ! (3*pi2*|rho|)^(1/3)
  ! |grho|
  ! |grho|/(2*kf*|rho|)
  ! n*ds/d(gn)
  ! exchange energy LDA part
  ! exchange energy gradient part
  real(DP) :: dfx, f1, f2
  ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
  real(DP) :: third, c1, c2, c5
  parameter (third = 1.0_dp / 3.0_dp, c1 = 0.75_dp / pi , &
       c2 = 3.093667726280136_dp, c5 = 4.0_dp * third)
  ! parameters of the functional
  real(DP) :: k (3), mu(3)
  data k / 0.804_dp, 1.245_dp, 0.804_dp /, &
       mu/ 0.21951_dp, 0.21951_dp, 0.12345679012345679012_dp  /
  !
  do i=1,length
     if ((rho(i).gt.small).and.(grho(i).gt.small**2)) then
        agrho = sqrt(grho(i))
        kf = c2 * rho(i)**third
        dsg = 0.5_dp / kf
        s1 = agrho * dsg / rho(i)
        !
        !   Energy
        f1 = s1*s1 * mu(iflag) / k(iflag)
        f2 = 1.0_dp / (1.0_dp + f1)
        fx = k(iflag) * (1.0_dp - f2)
        exunif = - c1 * kf
        sx(i) = exunif * fx
        !
        !   Potential
        dfx = 2.0_dp * mu(iflag) * s1 *f2*f2

        v1x(i) = sx(i) + exunif * (third * fx - c5 * dfx * s1)
        v2x(i) = exunif * dfx * dsg / agrho
        sx(i) = sx(i) * rho(i)

     else
        v1x(i) = 0.0_dp
        v2x(i) = 0.0_dp
        sx(i) = 0.0_dp
     end if
  end do

end subroutine pbex_vec
!
!---------------------------------------------------------------
subroutine pbec (rho, grho, iflag, sc, v1c, v2c)
  !---------------------------------------------------------------
  !
  ! PBE correlation (without LDA part)
  ! iflag=1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
  ! iflag=2: J.P.Perdew et al., PRL 100, 136406 (2008).
  ! iflag=3: L. Chiodo et al, PRL 108, 126402 (2012)  (PBEQ2D)
  !
  USE kinds, ONLY : DP
#ifdef __LIBXC
  use xc_f90_types_m
  use xc_f90_lib_m
#endif
  implicit none
  integer, intent(in) :: iflag
  real(DP), intent(in) :: rho, grho
  real(DP), intent(out):: sc, v1c, v2c
#ifdef __LIBXC
  ! local variables
  integer :: func_id = -1 ! not set
  integer :: size = 1
  TYPE(xc_f90_pointer_t) :: xc_func
  TYPE(xc_f90_pointer_t) :: xc_info
  real(dp) :: exc, ec_lda = 0.0d0 , vc_lda = 0.0d0

  if (iflag.eq.1)  func_id = 130
  if (iflag.eq.2)  func_id = 133
  if (iflag.eq.3)  call errore('pbec','case not implemented with libxc',iflag)

  call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)      
  call xc_f90_gga_exc_vxc(xc_func, size, rho, grho, exc, v1c, v2c)       
  call xc_f90_func_end(xc_func) 
  ! remove PW correlation for compatibility with QE  
  call xc_f90_func_init(xc_func, xc_info, 12, XC_UNPOLARIZED)    
  call xc_f90_lda_exc_vxc(xc_func,size , rho, ec_lda, vc_lda)
  call xc_f90_func_end(xc_func) 
  exc = exc - ec_lda 
  v1c = v1c - vc_lda
  
  sc = exc * rho  ! e_x = rho * \epsilon_x
  v2c = v2c*2.0_dp

#else
  real(DP), parameter :: ga = 0.031091d0
  real(DP) :: be (3)
!             pbe           pbesol   pbeq2d
  data be / 0.066725d0, 0.046d0,     0.066725d0/
  real(DP), parameter :: third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0
  real(DP), parameter :: xkf = 1.919158292677513d0, xks = 1.128379167095513d0
  ! pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
  real(DP) :: kf, ks, rs, ec, vc, t, expe, af, bf, y, xy, qy
  real(DP) :: s1, h0, dh0, ddh0, sc2D, v1c2D, v2c2D
  !
  rs = pi34 / rho**third
  call pw (rs, 1, ec, vc)
  kf = xkf / rs
  ks = xks * sqrt (kf)
  t = sqrt (grho) / (2.d0 * ks * rho)
  expe = exp ( - ec / ga)
  af = be(iflag) / ga * (1.d0 / (expe-1.d0) )
  bf = expe * (vc - ec)
  y = af * t * t
  xy = (1.d0 + y) / (1.d0 + y + y * y)
  qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
  s1 = 1.d0 + be(iflag) / ga * t * t * xy
  h0 = ga * log (s1)
  dh0 = be(iflag) * t * t / s1 * ( - 7.d0 / 3.d0 * xy - qy * (af * bf / &
       be(iflag)-7.d0 / 3.d0) )
  ddh0 = be(iflag) / (2.d0 * ks * ks * rho) * (xy - qy) / s1
  sc = rho * h0
  v1c = h0 + dh0
  v2c = ddh0
! q2D
  if (iflag == 3)then
     call cpbe2d(rho,grho,sc2D,v1c2D,v2c2D)
     sc=sc+sc2D
     v1c=v1c+v1c2D
     v2c=v2c+v2c2D
  endif
#endif
  !
  return
end subroutine pbec

!---------------------------------------------------------------
subroutine cpbe2d(rho,grho,sc,v1c,v2c)
!---------------------------------------------------------------

  ! 2D correction (last term of Eq. 5, PRL 108, 126402 (2012))

  USE kinds, ONLY : dp
  USE constants, ONLY : pi
  IMPLICIT NONE
  !
  REAL(dp), INTENT(in) :: rho, grho 
  REAL(dp), INTENT(OUT) :: sc, v1c, v2c
  !
  REAL(dp), PARAMETER:: ex1=0.333333333333333333_dp, ex2=1.166666666666667_dp
  REAL(dp), PARAMETER:: ex3=ex2+1.0_dp
  REAL(dp) :: fac1, fac2, zeta, phi, gr, rs, drsdn, akf, aks, t, dtdn, dtdgr
  REAL(dp) :: p, a, g, alpha1, beta1,beta2,beta3,beta4, dgdrs, epsc, depscdrs
  REAL(dp) :: c, gamma1, beta, aa, cg, adddepsc, h, dhdaa, dhdt, dhdrs
  REAL(dp) :: epscpbe, depscpbedrs, depscpbedt, a0,a1,a2, b0,b1,b2, c0,c1,c2
  REAL(dp) :: e0,e1,e2, f0,f1,f2, g0,g1,g2, h0,h1,h2, d0,d1,d2, ff, dffdt
  REAL(dp) :: rs3d, rs2d, drs2ddrs3d, eps2d, deps2ddrs2, depsGGAdrs, depsGGAdt
  REAL(dp) :: drs2ddt, rs2, ec, decdn, decdgr, daadepsc
  !
  fac1=(3.d0*pi*pi)**ex1
  fac2=sqrt(4.d0*fac1/pi)

  zeta=0.d0
  phi=1.d0
!
  gr=sqrt (grho)
!
  rs=(3.d0/4.d0/pi/rho)**ex1
  drsdn = -dble(3 ** (0.1D1 / 0.3D1)) * dble(2 ** (0.1D1 / 0.3D1)) * &
  0.3141592654D1 ** (-0.1D1 / 0.3D1) * (0.1D1 / rho) ** (-0.2D1 / &
  0.3D1) / rho ** 2 / 0.6D1
!
  akf=(3.d0*pi*pi*rho)**(1.d0/3.d0)
  aks=dsqrt(4.d0*akf/pi)
  t=gr/2.d0/phi/aks/rho
  dtdn=-7.d0/6.d0*gr/2.d0/phi/dsqrt(4.d0/pi)/   &
  ((3.d0*pi*pi)**(1.d0/6.d0))/(rho**(13.d0/6.d0))
  dtdgr=1.d0/2.d0/phi/aks/rho
!
! for the LDA correlation
  p=1.d0
  A=0.031091d0
  alpha1=0.21370d0
  beta1=7.5957d0
  beta2=3.5876d0
  beta3=1.6382d0
  beta4=0.49294d0
  G = -0.2D1 * A * dble(1 + alpha1 * rs) * log(0.1D1 + 0.1D1 / A / ( &
  beta1 * sqrt(dble(rs)) + dble(beta2 * rs) + dble(beta3 * rs ** ( &
  0.3D1 / 0.2D1)) + dble(beta4 * rs ** (p + 1))) / 0.2D1)

  dGdrs = -0.2D1 * A * alpha1 * log(0.1D1 + 0.1D1 / A / (beta1 * sqrt(rs) &
   + beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs ** &
  (p + 1)) / 0.2D1) + (0.1D1 + alpha1 * rs) / (beta1 * sqrt(rs) +  &
  beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1)) &
  ** 2 * (beta1 * rs ** (-0.1D1 / 0.2D1) / 0.2D1 + beta2 + 0.3D1 / &
  0.2D1 * beta3 * sqrt(rs) + beta4 * rs ** (p + 1) * dble(p + 1) /  &
  rs) / (0.1D1 + 0.1D1 / A / (beta1 * sqrt(rs) + beta2 * rs + beta3 * &
  rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1)) / 0.2D1)
!
  epsc=G
  depscdrs=dGdrs
!
! PBE
  c=1.d0
  gamma1=0.031091d0
  beta=0.066725d0
!
  AA = beta / gamma1 / (exp(-epsc / gamma1 / phi ** 3) - 0.1D1)
  cg = beta / gamma1 ** 2 / (exp(-epsc/ gamma1 / phi ** 3) - 0.1D1) &
  ** 2 / phi ** 3 * exp(-epsc / gamma1 / phi ** 3)
  dAAdepsc=cg
!
  if(t.le.10.d0)then
  H = dble(gamma1) * phi ** 3 * log(dble(1 + beta / gamma1 * t ** 2 &
  * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4)))
!
  dHdAA = gamma1 * phi ** 3 * (beta / gamma1 * t ** 4 / (1 + c * AA &
  * t ** 2 + AA ** 2 * t ** 4) - beta / gamma1 * t ** 2 * (1 + AA * &
  t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) ** 2 * (c * t **&
   2 + 2 * AA * t ** 4)) / (1 + beta / gamma1 * t ** 2 * (1 + AA * &
  t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4))
!
  dHdt = gamma1 * phi ** 3 * (2 * beta / gamma1 * t * (1 + AA * t ** &
   2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) + 2 * beta / gamma1 &
   * t ** 3 * AA / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) - beta / &
   gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** &
   2 * t ** 4) ** 2 * (2 * c * AA * t + 4 * AA ** 2 * t ** 3)) / (1 &
   + beta / gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t ** &
   2 + AA ** 2 * t ** 4))

  else

  H=gamma1*(phi**3)*dlog(1.d0+(beta/gamma1)*(1.d0/AA))
!
  dHdAA =gamma1*(phi**3)*1.d0/(1.d0+(beta/gamma1)*(1.d0/AA))* &
  (beta/gamma1)*(-1.d0/AA/AA)
!
  dHdt=0.d0
  endif
!
  dHdrs=dHdAA*dAAdepsc*depscdrs
!
  epscPBE=epsc+H
  depscPBEdrs=depscdrs+dHdrs
  depscPBEdt=dHdt
!
! START THE 2D CORRECTION 
!
  beta=1.3386d0
  a0=-0.1925d0
  a1=0.117331d0
  a2=0.0234188d0
  b0=0.0863136d0
  b1=-0.03394d0
  b2=-0.037093d0
  c0=0.057234d0
  c1=-0.00766765d0
  c2=0.0163618d0
  e0=1.0022d0
  e1=0.4133d0
  e2=1.424301d0
  f0=-0.02069d0
  f1=0.d0
  f2=0.d0
  g0=0.340d0
  g1=0.0668467d0
  g2=0.d0
  h0=0.01747d0
  h1=0.0007799d0
  h2=1.163099d0
  d0=-a0*h0
  d1=-a1*h1
  d2=-a2*h2
!
  ff = t ** 4 * (1 + t ** 2) / (1000000 + t ** 6)
  dffdt = 4 * t ** 3 * (1 + t ** 2) / (1000000 + t ** 6) + 2 * t ** &
  5 / (1000000 + t ** 6) - 6 * t ** 9 * (1 + t ** 2) / (1000000 + t &
  ** 6) ** 2
!
  rs3d=rs
  rs2d = 0.4552100000D0 * dble(3 ** (0.7D1 / 0.12D2)) * dble(4 ** ( &
  0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** ( &
  0.5D1 / 0.4D1) * sqrt(t)

  cg = 0.5690125000D0 * dble(3 ** (0.7D1 / 0.12D2)) * dble(4 ** ( &
  0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** (0.1D1 &
   / 0.4D1) * sqrt(t)
  drs2ddrs3d=cg

  cg = 0.2276050000D0 * dble(3 ** (0.7D1 / 0.12D2)) * dble(4 ** ( &
  0.5D1 / 0.12D2)) * dble((1 / pi) ** (-0.5D1 / 0.12D2)) * dble(rs3d ** &
   (0.5D1 / 0.4D1)) * dble(t ** (-0.1D1 / 0.2D1))
  drs2ddt=cg
  rs2=rs2d
!
  eps2d = (exp(-beta * rs2) - 0.1D1) * (-0.2D1 / 0.3D1 * sqrt(0.2D1) &
   * dble((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 /  &
  0.2D1)) / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * dble( &
  zeta ** 2) + 0.3D1 / 0.128D3 * dble(zeta ** 4)) * sqrt(0.2D1) / pi / &
   rs2) + a0 + (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 ** 3) * log(0.1D1 &
   + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** &
  2 + h0 * rs2 ** 3)) + (a1 + (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 ** &
   3) * log(0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1) &
   + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * dble(zeta ** 2) + (a2 + (b2  &
  * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) * log(0.1D1 + 0.1D1 / (e2 * &
   rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3 &
  ))) * dble(zeta ** 4)
!
  cg = -beta * exp(-beta * rs2) * (-0.2D1 / 0.3D1 * sqrt(0.2D1) *  &
  dble((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1)) &
  / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * dble(zeta ** &
   2) + 0.3D1 / 0.128D3 * dble(zeta ** 4)) * sqrt(0.2D1) / pi / rs2) &
   + (exp(-beta * rs2) - 0.1D1) * (0.2D1 / 0.3D1 * sqrt(0.2D1) * dble &
  ((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1)) / &
   pi / rs2 ** 2 - 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * dble(zeta  &
   ** 2) + 0.3D1 / 0.128D3 * dble(zeta ** 4)) * sqrt(0.2D1) / pi /  &
  rs2 ** 2) + (b0 + 0.2D1 * c0 * rs2 + 0.3D1 * d0 * rs2 ** 2) * log( &
  0.1D1 + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 &
   ** 2 + h0 * rs2 ** 3)) - (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 **  &
  3) / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0  &
  * rs2 ** 3) ** 2 * (e0 + 0.3D1 / 0.2D1 * f0 * sqrt(rs2) + 0.2D1 *  &
  g0 * rs2 + 0.3D1 * h0 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e0 * rs2 + f0 &
   * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0 * rs2 ** 3)) + (( &
  b1 + 0.2D1 * c1 * rs2 + 0.3D1 * d1 * rs2 ** 2) * log(0.1D1 + 0.1D1  &
  / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 *  &
  rs2 ** 3)) - (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 ** 3) / (e1 * rs2 &
   + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3) ** &
   2 * (e1 + 0.3D1 / 0.2D1 * f1 * sqrt(rs2) + 0.2D1 * g1 * rs2 +  &
  0.3D1 * h1 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** ( &
  0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * dble(zeta ** 2)  &
  + ((b2 + 0.2D1 * c2 * rs2 + 0.3D1 * d2 * rs2 ** 2) * log(0.1D1 +  &
  0.1D1 / (e2 * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 &
   * rs2 ** 3)) - (b2 * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) / (e2  &
  * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 **  &
  3) ** 2 * (e2 + 0.3D1 / 0.2D1 * f2 * sqrt(rs2) + 0.2D1 * g2 * rs2  &
  + 0.3D1 * h2 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e2 * rs2 + f2 * rs2 ** &
   (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3))) * dble(zeta ** &
   4)
  deps2ddrs2=cg
!
! GGA-2D
!
  depsGGAdrs=ff*(-depscPBEdrs+deps2ddrs2*drs2ddrs3d)
  depsGGAdt=dffdt*(-epscPBE+eps2d)+ff*  &
   (-depscPBEdt+deps2ddrs2*drs2ddt)
!
  ec=rho*(ff*(-epscPBE+eps2d))
!
  decdn=ff*(-epscPBE+eps2d)+rho*depsGGAdrs*drsdn+ &
      rho*depsGGAdt*dtdn
!
  decdgr=rho*depsGGAdt*dtdgr
!
  sc=ec
  v1c=decdn
  v2c=decdgr/gr
!
  RETURN
END subroutine cpbe2d
!
!---------------------------------------------------------------
subroutine sogga (rho, grho, sx, v1x, v2x)
  !-------------------------------------------------------------
  !
  ! SOGGA exchange
  ! 
  !
 USE kinds
 USE constants, ONLY : pi

  implicit none

  real(dp), intent(in)    :: rho, grho
  real(dp), intent(out)   :: sx, v1x, v2x
  ! input: charge and abs gradient
  ! output: energy
  ! output: potential

  ! local variables
  real(dp)                :: rho43, xs, xs2, dxs2_drho, dxs2_dgrho2
  real(dp)                :: CX, denom, C1, C2, Fso, Fpbe, ex, Fx, dFx_dxs2, dex_drho

  real(dp), parameter     :: one = 1.0_dp, two=2.0_dp, three = 3.0_dp,                &
  &                          four = 4.0_dp, eight = 8.0_dp,             &
  &                          f13 = one/three, f23 = two/three, f43 = four/three,      &
  &                          f34=three/four, f83 = eight/three, f12 = one/two

  real(dp), parameter     :: mu=0.12346d0, kapa=0.552d0
  !
  !_____________________________________________________________________


  CX    =  f34 * (three/pi)**f13            ! Cx LDA
  denom =  four * (three*pi**two)**f23
  C1    =  mu / denom
  C2    =  mu / (kapa * denom)

  rho43 = rho**f43
  xs    = grho / rho43
  xs2   = xs * xs

  dxs2_drho  = -f83 * xs2 / rho
  dxs2_dgrho2 = one /rho**f83

  ex        = - CX * rho43
  dex_drho  = - f43 * CX * rho**f13

  Fso       = kapa * (one - exp(-C2*xs2))
  Fpbe      = C1 * xs2 / (one + C2*xs2)

  Fx        =   f12 * (Fpbe + Fso)
  dFx_dxs2  =   f12 * (C1 / ((one + C2*xs2)**2) + C1*exp(-C2*xs2))

  !
  !   Energy
 !
  sx = Fx * ex
  !
  !   Potential
  !
  v1x = dex_drho * Fx  +  ex * dFx_dxs2 * dxs2_drho
  v2x = two * ex * dFx_dxs2 * dxs2_dgrho2

end subroutine sogga
!
!
!     ==================================================================
subroutine hcth(rho,grho,sx,v1x,v2x)
  !     ==================================================================
  !     HCTH/120, JCP 109, p. 6264 (1998)
  !     Parameters set-up after N.L. Doltsisnis & M. Sprik (1999)
  !     Present release: Mauro Boero, Tsukuba, 11/05/2004
  !--------------------------------------------------------------------------
  !     rhoa = rhob = 0.5 * rho
  !     grho is the SQUARE of the gradient of rho! --> gr=sqrt(grho)
  !     sx  : total exchange correlation energy at point r
  !     v1x : d(sx)/drho  (eq. dfdra = dfdrb in original)
  !     v2x : 1/gr*d(sx)/d(gr) (eq. 0.5 * dfdza = 0.5 * dfdzb in original)
  !--------------------------------------------------------------------------
  USE kinds, ONLY : DP
  USE constants, ONLY: pi
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x

  real(DP), parameter :: o3=1.0d0/3.0d0, o34=4.0d0/3.0d0, fr83=8.d0/3.d0
  real(DP) :: cg0(6), cg1(6), caa(6), cab(6), cx(6)
  real(DP) :: r3q2, r3pi, gr, rho_o3, rho_o34, xa, xa2, ra, rab, &
       dra_drho, drab_drho, g, dg, era1, dera1_dra, erab0, derab0_drab, &
       ex, dex_drho, uaa, uab, ux, ffaa, ffab,  dffaa_drho, dffab_drho,&
       denaa, denab, denx, f83rho, bygr, gaa, gab, gx, taa, tab, txx, &
       dgaa_drho, dgab_drho, dgx_drho, dgaa_dgr, dgab_dgr, dgx_dgr
  !
  r3q2=2.d0**(-o3)
  r3pi=(3.d0/pi)**o3
  !.....coefficients for pw correlation......................................
  cg0(1)= 0.031091d0
  cg0(2)= 0.213700d0
  cg0(3)= 7.595700d0
  cg0(4)= 3.587600d0
  cg0(5)= 1.638200d0
  cg0(6)= 0.492940d0
  cg1(1)= 0.015545d0
  cg1(2)= 0.205480d0
  cg1(3)=14.118900d0
  cg1(4)= 6.197700d0
  cg1(5)= 3.366200d0
  cg1(6)= 0.625170d0
  !......hcth-19-4.....................................
  caa(1)=  0.489508d+00
  caa(2)= -0.260699d+00
  caa(3)=  0.432917d+00
  caa(4)= -0.199247d+01
  caa(5)=  0.248531d+01
  caa(6)=  0.200000d+00
  cab(1)=  0.514730d+00
  cab(2)=  0.692982d+01
  cab(3)= -0.247073d+02
  cab(4)=  0.231098d+02
  cab(5)= -0.113234d+02
  cab(6)=  0.006000d+00
  cx(1) =  0.109163d+01
  cx(2) = -0.747215d+00
  cx(3) =  0.507833d+01
  cx(4) = -0.410746d+01
  cx(5) =  0.117173d+01
  cx(6)=   0.004000d+00
  !...........................................................................
  gr=DSQRT(grho)
  rho_o3=rho**(o3)
  rho_o34=rho**(o34)
  xa=1.25992105d0*gr/rho_o34
  xa2=xa*xa
  ra=0.781592642d0/rho_o3
  rab=r3q2*ra
  dra_drho=-0.260530881d0/rho_o34
  drab_drho=r3q2*dra_drho
  call pwcorr(ra,cg1,g,dg)
  era1=g
  dera1_dra=dg
  call pwcorr(rab,cg0,g,dg)
  erab0=g
  derab0_drab=dg
  ex=-0.75d0*r3pi*rho_o34
  dex_drho=-r3pi*rho_o3
  uaa=caa(6)*xa2
  uaa=uaa/(1.0d0+uaa)
  uab=cab(6)*xa2
  uab=uab/(1.0d0+uab)
  ux=cx(6)*xa2
  ux=ux/(1.0d0+ux)
  ffaa=rho*era1
  ffab=rho*erab0-ffaa
  dffaa_drho=era1+rho*dera1_dra*dra_drho
  dffab_drho=erab0+rho*derab0_drab*drab_drho-dffaa_drho
  ! mb-> i-loop removed
  denaa=1.d0/(1.0d0+caa(6)*xa2)
  denab=1.d0/(1.0d0+cab(6)*xa2)
  denx =1.d0/(1.0d0+cx(6)*xa2)
  f83rho=fr83/rho
  bygr=2.0d0/gr
  gaa=caa(1)+uaa*(caa(2)+uaa*(caa(3)+uaa*(caa(4)+uaa*caa(5))))
  gab=cab(1)+uab*(cab(2)+uab*(cab(3)+uab*(cab(4)+uab*cab(5))))
  gx=cx(1)+ux*(cx(2)+ux*(cx(3)+ux*(cx(4)+ux*cx(5))))
  taa=denaa*uaa*(caa(2)+uaa*(2.d0*caa(3)+uaa &
      *(3.d0*caa(4)+uaa*4.d0*caa(5))))
  tab=denab*uab*(cab(2)+uab*(2.d0*cab(3)+uab &
       *(3.d0*cab(4)+uab*4.d0*cab(5))))
  txx=denx*ux*(cx(2)+ux*(2.d0*cx(3)+ux &
       *(3.d0*cx(4)+ux*4.d0*cx(5))))
  dgaa_drho=-f83rho*taa
  dgab_drho=-f83rho*tab
  dgx_drho=-f83rho*txx
  dgaa_dgr=bygr*taa
  dgab_dgr=bygr*tab
  dgx_dgr=bygr*txx
  ! mb
  sx=ex*gx+ffaa*gaa+ffab*gab
  v1x=dex_drho*gx+ex*dgx_drho &
       +dffaa_drho*gaa+ffaa*dgaa_drho &
       +dffab_drho*gab+ffab*dgab_drho
  v2x=(ex*dgx_dgr+ffaa*dgaa_dgr+ffab*dgab_dgr)/gr
  return
end subroutine hcth
!-------------------------------------------------------------------=
subroutine pwcorr(r,c,g,dg)
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: r, g, dg, c(6)
  real(DP) :: r12, r32, r2, rb, drb, sb

  r12=dsqrt(r)
  r32=r*r12
  r2=r*r
  rb=c(3)*r12+c(4)*r+c(5)*r32+c(6)*r2
  sb=1.0d0+1.0d0/(2.0d0*c(1)*rb)
  g=-2.0d0*c(1)*(1.0d0+c(2)*r)*dlog(sb)
  drb=c(3)/(2.0d0*r12)+c(4)+1.5d0*c(5)*r12+2.0d0*c(6)*r
  dg=(1.0d0+c(2)*r)*drb/(rb*rb*sb)-2.0d0*c(1)*c(2)*dlog(sb)

  return
end subroutine pwcorr
!-----------------------------------------------------------------------------
!     ==================================================================
subroutine optx(rho,grho,sx,v1x,v2x)
!     OPTX, Handy et al. JCP 116, p. 5411 (2002) and refs. therein
!     Present release: Mauro Boero, Tsukuba, 10/9/2002
!--------------------------------------------------------------------------
!     rhoa = rhob = 0.5 * rho in LDA implementation
!     grho is the SQUARE of the gradient of rho! --> gr=sqrt(grho)
!     sx  : total exchange correlation energy at point r
!     v1x : d(sx)/drho
!     v2x : 1/gr*d(sx)/d(gr)
!--------------------------------------------------------------------------
  use kinds, only: DP
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x

  real(DP), parameter :: small=1.D-30, smal2=1.D-10
!.......coefficients and exponents....................
  real(DP), parameter :: o43=4.0d0/3.0d0, two13=1.259921049894873D0, &
       two53=3.174802103936399D0, gam=0.006D0, a1cx=0.9784571170284421D0,&
       a2=1.43169D0 
  real(DP) :: gr, rho43, xa, gamx2, uden, uu
  !.......OPTX in compact form..........................
  if(rho <= small) then
     sx=0.0D0
     v1x=0.0D0
     v2x=0.0D0
  else
     gr = max(grho,SMAL2)
     rho43=rho**o43
     xa=two13*DSQRT(gr)/rho43
     gamx2=gam*xa*xa
     uden=1.d+00/(1.d+00+gamx2)
     uu=a2*gamx2*gamx2*uden*uden
     uden=rho43*uu*uden
     sx=-rho43*(a1cx+uu)/two13
     v1x=o43*(sx+two53*uden)/rho
     v2x=-two53*uden/gr
  endif
  return
end subroutine optx
!
!---------------------------------------------------------------
subroutine wcx (rho, grho, sx, v1x, v2x)
  !---------------------------------------------------------------
  !
  ! Wu-Cohen exchange (without Slater exchange):
  ! Z. Wu and R. E. Cohen, PRB 73, 235116 (2006)
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  ! input: charge and squared gradient
  ! output: energy
  ! output: potential
  ! local variables
  real(DP) :: kf, agrho, s1, s2, es2, ds, dsg, exunif, fx
  ! (3*pi2*|rho|)^(1/3)
  ! |grho|
  ! |grho|/(2*kf*|rho|)
  ! s^2
  ! n*ds/dn
  ! n*ds/d(gn)
  ! exchange energy LDA part
  ! exchange energy gradient part
  real(DP) :: dxunif, dfx, f1, f2, f3, dfx1, x1, x2, x3, &
              dxds1, dxds2, dxds3
  ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) )
  real(DP) :: third, c1, c2, c5, c6, teneightyone
  parameter (third = 1.d0 / 3.d0, c1 = 0.75d0 / pi , &
       c2 = 3.093667726280136d0, c5 = 4.d0 * third, &
       teneightyone = 0.123456790123d0)
  ! parameters of the functional
  real(DP) :: k, mu, cwc
  parameter (k = 0.804d0, mu = 0.2195149727645171d0, cwc = 0.00793746933516d0) 
  !
  agrho = sqrt (grho)
  kf = c2 * rho**third
  dsg = 0.5d0 / kf
  s1 = agrho * dsg / rho
  s2 = s1 * s1
  es2 = exp(-s2)
  ds = - c5 * s1
  !
  !   Energy
  !
  ! x = 10/81 s^2 + (mu - 10/81) s^2 e^-s^2 + ln (1 + c s^4)
  x1 = teneightyone * s2 
  x2 = (mu - teneightyone) * s2 * es2
  x3 = log(1.d0 + cwc * s2 * s2)
  f1 = (x1 + x2 + x3) / k
  f2 = 1.d0 + f1
  f3 = k / f2
  fx = k - f3
  exunif = - c1 * kf
  sx = exunif * fx
  !
  !   Potential
  !
  dxunif = exunif * third
  dfx1 = f2 * f2
  dxds1 = teneightyone
  dxds2 = (mu - teneightyone) * es2 * (1.d0 - s2)
  dxds3 = 2.d0 * cwc * s2 / (1.d0 + cwc * s2 *s2)
  dfx = 2.d0 * s1 * (dxds1 + dxds2 + dxds3) / dfx1
  v1x = sx + dxunif * fx + exunif * dfx * ds
  v2x = exunif * dfx * dsg / agrho

  sx = sx * rho
  return
end subroutine wcx
!
!-----------------------------------------------------------------------
function dpz (rs, iflg)
  !-----------------------------------------------------------------------
  !  derivative of the correlation potential with respect to local density
  !  Perdew and Zunger parameterization of the Ceperley-Alder functional
  !
  use kinds, only: DP
  USE constants, ONLY: pi, fpi
  !
  implicit none
  !
  real(DP), intent (in) :: rs
  integer, intent(in) :: iflg
  real(DP) :: dpz
  !
  !  local variables
  !  a,b,c,d,gc,b1,b2 are the parameters defining the functional
  !
  real(DP), parameter :: a = 0.0311d0, b = -0.048d0, c = 0.0020d0, &
       d = -0.0116d0, gc = -0.1423d0, b1 = 1.0529d0, b2 = 0.3334d0,&
       a1 = 7.0d0 * b1 / 6.d0, a2 = 4.d0 * b2 / 3.d0
  real(DP) :: x, den, dmx, dmrs
  !
  !
  if (iflg == 1) then
     dmrs = a / rs + 2.d0 / 3.d0 * c * (log (rs) + 1.d0) + &
          (2.d0 * d-c) / 3.d0
  else
     x = sqrt (rs)
     den = 1.d0 + x * (b1 + x * b2)
     dmx = gc * ( (a1 + 2.d0 * a2 * x) * den - 2.d0 * (b1 + 2.d0 * &
          b2 * x) * (1.d0 + x * (a1 + x * a2) ) ) / den**3
     dmrs = 0.5d0 * dmx / x
  endif
  !
  dpz = - fpi * rs**4.d0 / 9.d0 * dmrs
  return
  !
end function dpz
!----------------------------------------------------------------------
!
! HSE (wPBE) stabbing starts HERE
!
! Note, that you can get PBEhole functional,
! M. Ernzerhof, J. Chem. Phys. 109, 3313 (1998),
! from this by just setting OMEGA=0
!
! These are wrappers to the reference implementation
!-----------------------------------------------------------------------
      SUBROUTINE pbexsr_lsd(RHOA,RHOB,GRHOAA,GRHOBB,sx, &
                            V1XA,V2XA,V1XB,V2XB,OMEGA)
!     ==--------------------------------------------------------------==
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(SMALL=1.D-20)
!     ==--------------------------------------------------------------==
      SXA=0.0D0
      SXB=0.0D0
      V1XA=0.0D0
      V2XA=0.0D0
      V1XB=0.0D0
      V2XB=0.0D0
      IF(RHOA.GT.SMALL.AND.GRHOAA.GT.SMALL) THEN
        CALL pbexsr(2.D0*RHOA, 4.D0*GRHOAA, SXA, V1XA, V2XA, OMEGA)
      ENDIF
      IF(RHOB.GT.SMALL.AND.GRHOBB.GT.SMALL) THEN
        CALL pbexsr(2.D0*RHOB, 4.D0*GRHOBB, SXB, V1XB, V2XB, OMEGA)
      ENDIF
      sx = 0.5D0*(SXA+SXB)
      V2XA = 2.D0*V2XA
      V2XB = 2.D0*V2XB          ! I HOPE THIS WORKS JUST LIKE THIS

!     ==--------------------------------------------------------------==
      RETURN
      END SUBROUTINE pbexsr_lsd
!
!-----------------------------------------------------------------------     
      SUBROUTINE pbexsr(RHO,GRHO,sx,V1X,V2X,OMEGA)
!-----------------------------------------------------------------------
!
!      INCLUDE 'cnst.inc'
      use kinds, ONLY : DP

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER(SMALL=1.D-20,SMAL2=1.D-08)
      PARAMETER(US=0.161620459673995492D0,AX=-0.738558766382022406D0, &
                UM=0.2195149727645171D0,UK=0.8040D0,UL=UM/UK)
      REAL(DP), PARAMETER :: f1 = -1.10783814957303361_DP, alpha = 2.0_DP/3.0_DP
!     ==--------------------------------------------------------------==

!      CALL XC(RHO,EX,EC,VX,VC)
      RS = RHO**(1.0_DP/3.0_DP)
      VX = (4.0_DP/3.0_DP)*f1*alpha*RS

!      AA    = DMAX1(GRHO,SMAL2)
      AA    = GRHO
!      RR    = RHO**(-4.0_DP/3.0_DP)
      RR    = 1.0_DP/(RHO*RS)
      EX    = AX/RR
      S2    = AA*RR*RR*US*US

      S = SQRT(S2)
      IF(S.GT.8.3D0) THEN
        S = 8.572844D0 - 18.796223D0/S2
      ENDIF
      CALL wpbe_analy_erfc_approx_grad(RHO,S,OMEGA,FX,D1X,D2X)
      sx = EX*FX        ! - EX
      DSDN = -4.D0/3.D0*S/RHO
      V1X = VX*FX + (DSDN*D2X+D1X)*EX   ! - VX
      DSDG = US*RR
      V2X = EX*1.D0/SQRT(AA)*DSDG*D2X

! NOTE, here sx is the total energy density,
! not just the gradient correction energy density as e.g. in pbex()
! And the same goes for the potentials V1X, V2X

!     ==--------------------------------------------------------------==
      RETURN
      END SUBROUTINE pbexsr
!
! gau-pbe in
!
!-----------------------------------------------------------------------
      SUBROUTINE pbexgau_lsd(RHOA,RHOB,GRHOAA,GRHOBB,sx, &
                            V1XA,V2XA,V1XB,V2XB,alpha_gau)
!     ==--------------------------------------------------------------==
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(SMALL=1.D-20)
!     ==--------------------------------------------------------------==
      SXA=0.0D0
      SXB=0.0D0
      V1XA=0.0D0
      V2XA=0.0D0
      V1XB=0.0D0
      V2XB=0.0D0
      IF(RHOA.GT.SMALL.AND.GRHOAA.GT.SMALL) THEN
        CALL pbexgau(2.D0*RHOA, 4.D0*GRHOAA, SXA, V1XA, V2XA, &
                                                   alpha_gau)
      ENDIF
      IF(RHOB.GT.SMALL.AND.GRHOBB.GT.SMALL) THEN
        CALL pbexgau(2.D0*RHOB, 4.D0*GRHOBB, SXB, V1XB, V2XB, &
                                                   alpha_gau)
      ENDIF
      sx = 0.5D0*(SXA+SXB)
      V2XA = 2.D0*V2XA
      V2XB = 2.D0*V2XB          ! I HOPE THIS WORKS JUST LIKE THIS

!     ==--------------------------------------------------------------==
      RETURN
      END SUBROUTINE pbexgau_lsd
!
!-----------------------------------------------------------------------
      SUBROUTINE pbexgau(RHO,GRHO,sxsr,v1xsr,v2xsr,alpha_gau)
!-----------------------------------------------------------------------
!
      use kinds, ONLY : DP

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER(SMALL=1.D-20,SMAL2=1.D-08)
      PARAMETER(US=0.161620459673995492D0,AX=-0.738558766382022406D0, &
                UM=0.2195149727645171D0,UK=0.8040D0,UL=UM/UK)
      REAL(DP), PARAMETER :: f1 = -1.10783814957303361_DP, alpha = 2.0_DP/3.0_DP
!     ==--------------------------------------------------------------==

      RS = RHO**(1.0_DP/3.0_DP)
      VX = (4.0_DP/3.0_DP)*f1*alpha*RS
      AA    = GRHO
      RR    = 1.0_DP/(RHO*RS)
      EX    = AX/RR
! AX is 3/4/PI*(3*PI*PI)**(1/3). This is the same as -c1*c2 in pbex().
      S2    = AA*RR*RR*US*US
      S = SQRT(S2)
      IF(S.GT.10.D0) THEN
        S = 10.D0
      ENDIF
      CALL pbe_gauscheme(RHO,S,alpha_gau,FX,D1X,D2X)
      sxsr = EX*FX        ! - EX
      DSDN = -4.D0/3.D0*S/RHO
      V1Xsr = VX*FX + (DSDN*D2X+D1X)*EX   ! - VX
      DSDG = US*RR
      V2Xsr = EX*1.D0/SQRT(AA)*DSDG*D2X

! NOTE, here sx is the total energy density,
! not just the gradient correction energy density as e.g. in pbex()
! And the same goes for the potentials V1X, V2X

!     ==--------------------------------------------------------------==
      RETURN
      END SUBROUTINE pbexgau
!
!-----------------------------------------------------------------------
      SUBROUTINE pbe_gauscheme(rho,s,alpha_gau,Fx,dFxdr,dFxds)
!--------------------------------------------------------------------

      Implicit None
      Real*8 rho,s,alpha_gau,Fx,dFxdr,dFxds
!     input: charge and squared gradient and alpha_gau
!     output: GGA enhancement factor of gau-PBE
!     output: d(Fx)/d(s) , d(Fx)/d(rho)

      Real*8 Kx, Nx
!     PBE96 GGA enhancement factor
!     GGA enhancement factor of Gaussian Function

      Real*8 bx, cx, PI, sqrtpial, Prefac, term_PBE, Third, KsF
      Real*8 d1sdr, d1Kxds, d1Kxdr, d1bxdr, d1bxds, d1bxdKx, &
           d1Nxdbx,d1Nxdr, d1Nxds

      Real*8, external :: qe_erf,TayExp

      Real*8 Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten

      Save Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten
      Data Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten &
        / 0D0,1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0,10D0 /

      Real*8 k , mu
      Data k / 0.804d0 / , mu / 0.21951d0 /
!     parameters of PBE functional

      Third = One/Three
      PI = ACos(-One)
      KsF = (Three*PI*PI*rho)**Third
      sqrtpial = sqrt(PI/alpha_gau)
      Prefac = Two *sqrt(PI/alpha_gau) / Three

!     PBE96 GGA enhancement factor part
      term_PBE = One / (One + s*s*mu/k)
      Kx =  One + k - k * term_PBE

!     GGA enhancement factor of Gaussian Function part
      bx = sqrt(Kx*alpha_gau) / KsF

!      cx = exp(-One/Four/bx/bx) - One
      If(Abs(One/bx/bx) .lt. 1.0D-4) then
         cx = TayExp(-One/bx/bx)
      else
         cx = exp(-One/bx/bx) - One
      endIf

      Nx = bx * Prefac * ( sqrt(PI) * qe_erf(One/bx) + &
       (bx - Two*bx*bx*bx)*cx - Two*bx )

! for convergency
      If(Abs(Nx) .lt. 1.0D-15)then
        Nx = Zero
      else if ((One - Abs(Nx)) .lt. 1.0D-15)then
        Nx = One
      else
        Nx = Nx
      endIf
! for convergency end

      Fx =  Kx * Nx

!     1st derivatives
      d1sdr = - Four / Three * s / rho

      d1Kxds = Two * s * mu * term_PBE * term_PBE
      d1Kxdr = d1Kxds * d1sdr
      d1bxdKx = bx / (Two* Kx)

      d1bxdr = - bx /(Three*rho) + d1Kxdr * d1bxdKx

      d1bxds =  d1bxdKx * d1Kxds

      d1Nxdbx =  Nx/bx - Prefac * bx * Three * &
                  ( cx*(One + Two*bx*bx) + Two )

      d1Nxdr = d1Nxdbx * d1bxdr
      d1Nxds = d1Nxdbx * d1bxds

      dFxdr = d1Kxdr * Nx + Kx * d1Nxdr
      dFxds = d1Kxds * Nx + Kx * d1Nxds

      RETURN
      END SUBROUTINE pbe_gauscheme
!
      FUNCTION TayExp(X)
      Real*8 TAYEXP,X
      INTEGER NTERM,I
      Real*8 SUMVAL,IVAL,COEF
      PARAMETER (NTERM=16)

      SUMVAL = X
      IVAL = X
      COEF = 1.0D0
      DO 10 I = 2,NTERM
         COEF = COEF * I
         IVAL = IVAL * (X / COEF)
         SUMVAL = SUMVAL + IVAL
 10         CONTINUE
      TAYEXP = SUMVAL
      RETURN
      END FUNCTION TayExp
! gau-pbe out
!
!-----------------------------------------------------------------------
      SUBROUTINE wpbe_analy_erfc_approx_grad(rho,s,omega,Fx_wpbe, &
                      d1rfx,d1sfx)
!--------------------------------------------------------------------
!
!     wPBE Enhancement Factor (erfc approx.,analytical, gradients)
!
!--------------------------------------------------------------------

      Implicit None

      Real*8 rho,s,omega,Fx_wpbe,d1sfx,d1rfx

      Real*8 f12,f13,f14,f18,f23,f43,f32,f72,f34,f94,f1516,f98
      Real*8 pi,pi2,pi_23,srpi
      Real*8 Three_13

      Real*8 ea1,ea2,ea3,ea4,ea5,ea6,ea7,ea8
      Real*8 eb1
      Real*8 A,B,C,D,E
      Real*8 Ha1,Ha2,Ha3,Ha4,Ha5
      Real*8 Fc1,Fc2
      Real*8 EGa1,EGa2,EGa3
      Real*8 EGscut,wcutoff,expfcutoff

      Real*8 xkf, xkfrho
      Real*8 w,w2,w3,w4,w5,w6,w7,w8
      Real*8 d1rw
      Real*8 A2,A3,A4,A12,A32,A52,A72
      Real*8 X
      Real*8 s2,s3,s4,s5,s6

      Real*8 H,F
      Real*8 Hnum,Hden,d1sHnum,d1sHden
      Real*8 d1sH,d1sF
      Real*8 G_a,G_b,EG
      Real*8 d1sG_a,d1sG_b,d1sEG

      Real*8 Hsbw,Hsbw2,Hsbw3,Hsbw4,Hsbw12,Hsbw32,Hsbw52,Hsbw72
      Real*8 DHsbw,DHsbw2,DHsbw3,DHsbw4,DHsbw5
      Real*8 DHsbw12,DHsbw32,DHsbw52,DHsbw72,DHsbw92
      Real*8 d1sHsbw,d1rHsbw
      Real*8 d1sDHsbw,d1rDHsbw
      Real*8 HsbwA94,HsbwA9412
      Real*8 HsbwA942,HsbwA943,HsbwA945
      Real*8 piexperf,expei
      Real*8 piexperfd1,expeid1
      Real*8 d1spiexperf,d1sexpei
      Real*8 d1rpiexperf,d1rexpei
      Real*8 expei1,expei2,expei3,expei4

      Real*8 DHs,DHs2,DHs3,DHs4,DHs72,DHs92,DHsw,DHsw2,DHsw52,DHsw72
      Real*8 d1sDHs,d1rDHsw

      Real*8 np1,np2
      Real*8 d1rnp1,d1rnp2
      Real*8 t1,t2t9,t10,t10d1
      Real*8 f2,f3,f4,f5,f6,f7,f8,f9
      Real*8 f2d1,f3d1,f4d1,f5d1,f6d1,f8d1,f9d1
      Real*8 d1sf2,d1sf3,d1sf4,d1sf5,d1sf6,d1sf7,d1sf8,d1sf9
      Real*8 d1rf2,d1rf3,d1rf4,d1rf5,d1rf6,d1rf7,d1rf8,d1rf9
      Real*8 d1st1,d1rt1
      Real*8 d1st2t9,d1rt2t9
      Real*8 d1st10,d1rt10
      Real*8 d1sterm1,d1rterm1,term1d1
      Real*8 d1sterm2
      Real*8 d1sterm3,d1rterm3
      Real*8 d1sterm4,d1rterm4
      Real*8 d1sterm5,d1rterm5

      Real*8 term1,term2,term3,term4,term5

      Real*8 ax,um,uk,ul
      Real*8 gc1,gc2

      Real*8, external :: qe_erf, qe_erfc
!      Real*8 ei
      Real*8, external :: expint

      Real*8 Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten
      Real*8 Fifteen,Sixteen
      Real*8 r12,r64,r36,r81,r256,r384,r864,r1944,r4374
      Real*8 r20,r25,r27,r48,r120,r128,r144,r288,r324,r512,r729
      Real*8 r30,r32,r75,r243,r2187,r6561,r40,r105,r54,r135
      Real*8 r1215,r15309

      Save Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten
      Data Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten &
        / 0D0,1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0,10D0 /
      Save Fifteen,Sixteen
      Data Fifteen,Sixteen / 1.5D1, 1.6D1 /
      Save r36,r64,r81,r256,r384,r864,r1944,r4374
      Data r36,r64,r81,r256,r384,r864,r1944,r4374 &
        / 3.6D1,6.4D1,8.1D1,2.56D2,3.84D2,8.64D2,1.944D3,4.374D3 /
      Save r27,r48,r120,r128,r144,r288,r324,r512,r729
      Data r27,r48,r120,r128,r144,r288,r324,r512,r729 &
        / 2.7D1,4.8D1,1.2D2,1.28D2,1.44D2,2.88D2,3.24D2,5.12D2,7.29D2 /
      Save r20,r32,r243,r2187,r6561,r40
      Data r20,r32,r243,r2187,r6561,r40 &
        / 2.0d1,3.2D1,2.43D2,2.187D3,6.561D3,4.0d1 /
      Save r12,r25,r30,r54,r75,r105,r135,r1215,r15309
      Data r12,r25,r30,r54,r75,r105,r135,r1215,r15309 &
        / 1.2D1,2.5d1,3.0d1,5.4D1,7.5d1,1.05D2,1.35D2,1.215D3,1.5309D4 /

!     General constants

      f12    = 0.5d0
      f13    = One/Three
      f14    = 0.25d0
      f18    = 0.125d0

      f23    = Two * f13
      f43    = Two * f23

      f32    = 1.5d0
      f72    = 3.5d0
      f34    = 0.75d0
      f94    = 2.25d0
      f98    = 1.125d0
      f1516  = Fifteen / Sixteen

      pi     = ACos(-One)
      pi2    = pi*pi
      pi_23  = pi2**f13
      srpi   = sqrt(pi)

      Three_13 = Three**f13

!     Constants from fit

      ea1 = -1.128223946706117d0
      ea2 = 1.452736265762971d0
      ea3 = -1.243162299390327d0
      ea4 = 0.971824836115601d0
      ea5 = -0.568861079687373d0
      ea6 = 0.246880514820192d0
      ea7 = -0.065032363850763d0
      ea8 = 0.008401793031216d0

      eb1 = 1.455915450052607d0

!     Constants for PBE hole

      A      =  1.0161144d0
      B      = -3.7170836d-1
      C      = -7.7215461d-2
      D      =  5.7786348d-1
      E      = -5.1955731d-2
      X      = - Eight/Nine

!     Constants for fit of H(s) (PBE)

      Ha1    = 9.79681d-3
      Ha2    = 4.10834d-2
      Ha3    = 1.87440d-1
      Ha4    = 1.20824d-3
      Ha5    = 3.47188d-2

!     Constants for F(H) (PBE)

      Fc1    = 6.4753871d0
      Fc2    = 4.7965830d-1

!     Constants for polynomial expansion for EG for small s

      EGa1   = -2.628417880d-2
      EGa2   = -7.117647788d-2
      EGa3   =  8.534541323d-2

!     Constants for large x expansion of exp(x)*ei(-x)

      expei1 = 4.03640D0
      expei2 = 1.15198D0
      expei3 = 5.03627D0
      expei4 = 4.19160D0

!     Cutoff criterion below which to use polynomial expansion

      EGscut     = 8.0d-2
      wcutoff    = 1.4D1
      expfcutoff = 7.0D2

!     Calculate prelim variables

      xkf    = (Three*pi2*rho) ** f13
      xkfrho = xkf * rho

      A2 = A*A
      A3 = A2*A
      A4 = A3*A
      A12 = Sqrt(A)
      A32 = A12*A
      A52 = A32*A
      A72 = A52*A

      w      = omega / xkf
      w2    = w * w
      w3    = w2 * w
      w4    = w2 * w2
      w5    = w3 * w2
      w6    = w5 * w
      w7    = w6 * w
      w8    = w7 * w

      d1rw  = -(One/(Three*rho))*w

      X      = - Eight/Nine

      s2     = s*s
      s3     = s2*s
      s4     = s2*s2
      s5     = s4*s
      s6     = s5*s

!     Calculate wPBE enhancement factor

      Hnum    = Ha1*s2 + Ha2*s4
      Hden    = One + Ha3*s4 + Ha4*s5 + Ha5*s6

      H       = Hnum/Hden

      d1sHnum = Two*Ha1*s + Four*Ha2*s3
      d1sHden = Four*Ha3*s3 + Five*Ha4*s4 + Six*Ha5*s5

      d1sH    = (Hden*d1sHnum - Hnum*d1sHden) / (Hden*Hden)

      F      = Fc1*H + Fc2
      d1sF   = Fc1*d1sH

!     Change exponent of Gaussian if we're using the simple approx.

      if(w .gt. wcutoff) then

        eb1 = 2.0d0

      endif

!     Calculate helper variables (should be moved later on...)

      Hsbw = s2*H + eb1*w2
      Hsbw2 = Hsbw*Hsbw
      Hsbw3 = Hsbw2*Hsbw
      Hsbw4 = Hsbw3*Hsbw
      Hsbw12 = Sqrt(Hsbw)
      Hsbw32 = Hsbw12*Hsbw
      Hsbw52 = Hsbw32*Hsbw
      Hsbw72 = Hsbw52*Hsbw

      d1sHsbw  = d1sH*s2 + Two*s*H
      d1rHsbw  = Two*eb1*d1rw*w

      DHsbw = D + s2*H + eb1*w2
      DHsbw2 = DHsbw*DHsbw
      DHsbw3 = DHsbw2*DHsbw
      DHsbw4 = DHsbw3*DHsbw
      DHsbw5 = DHsbw4*DHsbw
      DHsbw12 = Sqrt(DHsbw)
      DHsbw32 = DHsbw12*DHsbw
      DHsbw52 = DHsbw32*DHsbw
      DHsbw72 = DHsbw52*DHsbw
      DHsbw92 = DHsbw72*DHsbw

      HsbwA94   = f94 * Hsbw / A
      HsbwA942  = HsbwA94*HsbwA94
      HsbwA943  = HsbwA942*HsbwA94
      HsbwA945  = HsbwA943*HsbwA942
      HsbwA9412 = Sqrt(HsbwA94)

      DHs    = D + s2*H
      DHs2   = DHs*DHs
      DHs3   = DHs2*DHs
      DHs4   = DHs3*DHs
      DHs72  = DHs3*sqrt(DHs)
      DHs92  = DHs72*DHs

      d1sDHs = Two*s*H + s2*d1sH

      DHsw   = DHs + w2
      DHsw2  = DHsw*DHsw
      DHsw52 = sqrt(DHsw)*DHsw2
      DHsw72 = DHsw52*DHsw

      d1rDHsw = Two*d1rw*w

      if(s .gt. EGscut) then

        G_a    = srpi * (Fifteen*E + Six*C*(One+F*s2)*DHs + &
                         Four*B*(DHs2) + Eight*A*(DHs3))    &
                      * (One / (Sixteen * DHs72))           &
                       - f34*pi*sqrt(A) * exp(f94*H*s2/A) * &
                         (One - qe_erf(f32*s*sqrt(H/A)))

        d1sG_a = (One/r32)*srpi *                           &
                 ((r36*(Two*H + d1sH*s) / (A12*sqrt(H/A)))  &
                  + (One/DHs92) *                           &
                     (-Eight*A*d1sDHs*DHs3 - r105*d1sDHs*E  &
                      -r30*C*d1sDHs*DHs*(One+s2*F)          &
                      +r12*DHs2*(-B*d1sDHs + C*s*(d1sF*s + Two*F)))  &
                  - ((r54*exp(f94*H*s2/A)*srpi*s*(Two*H+d1sH*s)*     &
                     qe_erfc(f32*sqrt(H/A)*s))                         &
                     / A12))

        G_b    = (f1516 * srpi * s2) / DHs72

        d1sG_b = (Fifteen*srpi*s*(Four*DHs - Seven*d1sDHs*s)) &
                 / (r32*DHs92)

        EG     = - (f34*pi + G_a) / G_b

        d1sEG  = (-Four*d1sG_a*G_b + d1sG_b*(Four*G_a + Three*pi)) &
                 / (Four*G_b*G_b)

      else

        EG    = EGa1 + EGa2*s2 + EGa3*s4
        d1sEG = Two*EGa2*s + Four*EGa3*s3

      endif

!    Calculate the terms needed in any case


      term2 =       (DHs2*B + DHs*C + Two*E + DHs*s2*C*F + Two*s2*EG) / &
                    (Two*DHs3)

      d1sterm2 = (-Six*d1sDHs*(EG*s2 + E)                     &
                  + DHs2 * (-d1sDHs*B + s*C*(d1sF*s + Two*F)) &
                  + Two*DHs * (Two*EG*s - d1sDHs*C            &
                  + s2 * (d1sEG - d1sDHs*C*F)))               &
                 / (Two*DHs4)

      term3 = - w  * (Four*DHsw2*B + Six*DHsw*C + Fifteen*E &
                      + Six*DHsw*s2*C*F + Fifteen*s2*EG) /  &
                     (Eight*DHs*DHsw52)

      d1sterm3 = w * (Two*d1sDHs*DHsw * (Four*DHsw2*B         &
                         + Six*DHsw*C + Fifteen*E             &
                         + Three*s2*(Five*EG + Two*DHsw*C*F)) &
                      + DHs * (r75*d1sDHs*(EG*s2 + E)         &
                         + Four*DHsw2*(d1sDHs*B               &
                              - Three*s*C*(d1sF*s + Two*F))   &
                         - Six*DHsw*(-Three*d1sDHs*C          &
                              + s*(Ten*EG + Five*d1sEG*s      &
                                  - Three*d1sDHs*s*C*F))))    &
                 / (Sixteen*DHs2*DHsw72)

      d1rterm3 = (-Two*d1rw*DHsw * (Four*DHsw2*B              &
                         + Six*DHsw*C + Fifteen*E             &
                         + Three*s2*(Five*EG + Two*DHsw*C*F)) &
                      + w * d1rDHsw * (r75*(EG*s2 + E)        &
                         + Two*DHsw*(Two*DHsw*B + Nine*C      &
                                     + Nine*s2*C*F)))         &
                 / (Sixteen*DHs*DHsw72)

      term4 = - w3 * (DHsw*C + Five*E + DHsw*s2*C*F + Five*s2*EG) /  &
                     (Two*DHs2*DHsw52)

      d1sterm4 = (w3 * (Four*d1sDHs*DHsw * (DHsw*C + Five*E   &
                             + s2 * (Five*EG + DHsw*C*F))     &
                        + DHs * (r25*d1sDHs*(EG*s2 + E)       &
                             - Two*DHsw2*s*C*(d1sF*s + Two*F) &
                             + DHsw * (Three*d1sDHs*C + s*(-r20*EG  &
                                   - Ten*d1sEG*s              &
                                   + Three*d1sDHs*s*C*F)))))  &
                 / (Four*DHs3*DHsw72)

      d1rterm4 = (w2 * (-Six*d1rw*DHsw * (DHsw*C + Five*E   &
                             + s2 * (Five*EG + DHsw*C*F))   &
                        + w * d1rDHsw * (r25*(EG*s2 + E) +  &
                             Three*DHsw*C*(One + s2*F))))  &
                 / (Four*DHs2*DHsw72)

      term5 = - w5 * (E + s2*EG) / &
                     (DHs3*DHsw52)

      d1sterm5 = (w5 * (Six*d1sDHs*DHsw*(EG*s2 + E)               &
                        + DHs * (-Two*DHsw*s * (Two*EG + d1sEG*s) &
                             + Five*d1sDHs * (EG*s2 + E))))       &
                 / (Two*DHs4*DHsw72)

      d1rterm5 = (w4 * Five*(EG*s2 + E) * (-Two*d1rw*DHsw   &
                                           + d1rDHsw * w))  &
                 / (Two*DHs3*DHsw72)


      if((s.gt.0.0d0).or.(w.gt.0.0d0)) then

        t10    = (f12)*A*Log(Hsbw / DHsbw)
        t10d1  = f12*A*(One/Hsbw - One/DHsbw)
        d1st10 = d1sHsbw*t10d1
        d1rt10 = d1rHsbw*t10d1

      endif

!     Calculate exp(x)*f(x) depending on size of x

      if(HsbwA94 .lt. expfcutoff) then

        piexperf = pi*Exp(HsbwA94)*qe_erfc(HsbwA9412)
!        expei    = Exp(HsbwA94)*Ei(-HsbwA94)
        expei    = Exp(HsbwA94)*(-expint(1,HsbwA94))

      else

!       print *,rho,s," LARGE HsbwA94"

        piexperf = pi*(One/(srpi*HsbwA9412)          &
                   - One/(Two*Sqrt(pi*HsbwA943))     &
                   + Three/(Four*Sqrt(pi*HsbwA945)))

        expei  = - (One/HsbwA94) *                         &
                   (HsbwA942 + expei1*HsbwA94 + expei2) /  &
                   (HsbwA942 + expei3*HsbwA94 + expei4)

      endif

!     Calculate the derivatives (based on the orig. expression)
!     --> Is this ok? ==> seems to be ok...

      piexperfd1  = - (Three*srpi*sqrt(Hsbw/A))/(Two*Hsbw)  &
                    + (Nine*piexperf)/(Four*A)
      d1spiexperf = d1sHsbw*piexperfd1
      d1rpiexperf = d1rHsbw*piexperfd1

      expeid1  = f14*(Four/Hsbw + (Nine*expei)/A)
      d1sexpei = d1sHsbw*expeid1
      d1rexpei = d1rHsbw*expeid1

      if (w .eq. Zero) then

!       Fall back to original expression for the PBE hole

        t1 = -f12*A*expei
        d1st1 = -f12*A*d1sexpei
        d1rt1 = -f12*A*d1rexpei

!        write(*,*) s, t1, t10, d1st1,d1rt1,d1rt10

        if(s .gt. 0.0D0) then

          term1    = t1 + t10
          d1sterm1 = d1st1 + d1st10
          d1rterm1 = d1rt1 + d1rt10

          Fx_wpbe = X * (term1 + term2)

          d1sfx = X * (d1sterm1 + d1sterm2)
          d1rfx = X * d1rterm1

        else

          Fx_wpbe = 1.0d0

! TODO    This is checked to be true for term1
!         How about the other terms???

          d1sfx   = 0.0d0
          d1rfx   = 0.0d0

        endif


      elseif(w .gt. wcutoff) then

!       Use simple Gaussian approximation for large w

!       print *,rho,s," LARGE w"

        term1 = -f12*A*(expei+log(DHsbw)-log(Hsbw))

        term1d1  = - A/(Two*DHsbw) - f98*expei
        d1sterm1 = d1sHsbw*term1d1
        d1rterm1 = d1rHsbw*term1d1

        Fx_wpbe = X * (term1 + term2 + term3 + term4 + term5)

        d1sfx = X * (d1sterm1 + d1sterm2 + d1sterm3  &
                              + d1sterm4 + d1sterm5)

        d1rfx = X * (d1rterm1 + d1rterm3 + d1rterm4 + d1rterm5)

      else

!       For everything else, use the full blown expression

!       First, we calculate the polynomials for the first term

         np1    = -f32*ea1*A12*w + r27*ea3*w3/(Eight*A12)     &
                  - r243*ea5*w5/(r32*A32) + r2187*ea7*w7/(r128*A52)

        d1rnp1 = - f32*ea1*d1rw*A12 + (r81*ea3*d1rw*w2)/(Eight*A12) &
                 - (r1215*ea5*d1rw*w4)/(r32*A32)                    &
                 + (r15309*ea7*d1rw*w6)/(r128*A52)

        np2 = -A + f94*ea2*w2 - r81*ea4*w4/(Sixteen*A)        &
              + r729*ea6*w6/(r64*A2) - r6561*ea8*w8/(r256*A3)


        d1rnp2 =   f12*(Nine*ea2*d1rw*w)         &
                 - (r81*ea4*d1rw*w3)/(Four*A)    &
                 + (r2187*ea6*d1rw*w5)/(r32*A2)  &
                 - (r6561*ea8*d1rw*w7)/(r32*A3)

!       The first term is

        t1    = f12*(np1*piexperf + np2*expei)
        d1st1 = f12*(d1spiexperf*np1 + d1sexpei*np2)
        d1rt1 = f12*(d1rnp2*expei + d1rpiexperf*np1 +  &
                     d1rexpei*np2 + d1rnp1*piexperf)

!       The factors for the main polynomoal in w and their derivatives

        f2    = (f12)*ea1*srpi*A / DHsbw12
        f2d1  = - ea1*srpi*A / (Four*DHsbw32)
        d1sf2 = d1sHsbw*f2d1
        d1rf2 = d1rHsbw*f2d1

        f3    = (f12)*ea2*A / DHsbw
        f3d1  = - ea2*A / (Two*DHsbw2)
        d1sf3 = d1sHsbw*f3d1
        d1rf3 = d1rHsbw*f3d1

        f4    =  ea3*srpi*(-f98 / Hsbw12     &
                 + f14*A / DHsbw32)
        f4d1  = ea3*srpi*((Nine/(Sixteen*Hsbw32))-   &
                          (Three*A/(Eight*DHsbw52)))
        d1sf4 = d1sHsbw*f4d1
        d1rf4 = d1rHsbw*f4d1

        f5    = ea4*(One/r128) * (-r144*(One/Hsbw)   &
                 + r64*(One/DHsbw2)*A)
        f5d1  = ea4*((f98/Hsbw2)-(A/DHsbw3))
        d1sf5 = d1sHsbw*f5d1
        d1rf5 = d1rHsbw*f5d1

        f6    = ea5*(Three*srpi*(Three*DHsbw52*(Nine*Hsbw-Two*A) &
                 + Four*Hsbw32*A2))                              &
                 / (r32*DHsbw52*Hsbw32*A)
        f6d1  = ea5*srpi*((r27/(r32*Hsbw52))-        &
                    (r81/(r64*Hsbw32*A))-            &
                    ((Fifteen*A)/(Sixteen*DHsbw72)))
        d1sf6 = d1sHsbw*f6d1
        d1rf6 = d1rHsbw*f6d1

        f7    = ea6*(((r32*A)/DHsbw3                 &
                 + (-r36 + (r81*s2*H)/A)/Hsbw2)) / r32
        d1sf7 = ea6*(Three*(r27*d1sH*DHsbw4*Hsbw*s2 +           &
                Eight*d1sHsbw*A*(Three*DHsbw4 - Four*Hsbw3*A) + &
                r54*DHsbw4*s*(Hsbw - d1sHsbw*s)*H))/            &
                (r32*DHsbw4*Hsbw3*A)
        d1rf7 = ea6*d1rHsbw*((f94/Hsbw3)-((Three*A)/DHsbw4)     &
                           -((r81*s2*H)/(Sixteen*Hsbw3*A)))

        f8    = ea7*(-Three*srpi*(-r40*Hsbw52*A3                &
                 +Nine*DHsbw72*(r27*Hsbw2-Six*Hsbw*A+Four*A2))) &
                 / (r128 * DHsbw72*Hsbw52*A2)
        f8d1  = ea7*srpi*((r135/(r64*Hsbw72)) + (r729/(r256*Hsbw32*A2))  &
                         -(r243/(r128*Hsbw52*A))                         &
                         -((r105*A)/(r32*DHsbw92)))
        d1sf8 = d1sHsbw*f8d1
        d1rf8 = d1rHsbw*f8d1

        f9    = (r324*ea6*eb1*DHsbw4*Hsbw*A                      &
                + ea8*(r384*Hsbw3*A3 + DHsbw4*(-r729*Hsbw2       &
                + r324*Hsbw*A - r288*A2))) / (r128*DHsbw4*Hsbw3*A2)
        f9d1  = -((r81*ea6*eb1)/(Sixteen*Hsbw3*A))               &
                + ea8*((r27/(Four*Hsbw4))+(r729/(r128*Hsbw2*A2)) &
                      -(r81/(Sixteen*Hsbw3*A))                   &
                      -((r12*A/DHsbw5)))
        d1sf9 = d1sHsbw*f9d1
        d1rf9 = d1rHsbw*f9d1

        t2t9    = f2*w  + f3*w2 + f4*w3 + f5*w4 + f6*w5          &
                        + f7*w6 + f8*w7 + f9*w8
        d1st2t9 = d1sf2*w + d1sf3*w2 + d1sf4*w3 + d1sf5*w4       &
                          + d1sf6*w5 + d1sf7*w6 + d1sf8*w7       &
                          + d1sf9*w8
        d1rt2t9 = d1rw*f2 + d1rf2*w + Two*d1rw*f3*w   &
                  + d1rf3*w2 + Three*d1rw*f4*w2       &
                  + d1rf4*w3 + Four*d1rw*f5*w3        &
                  + d1rf5*w4 + Five*d1rw*f6*w4        &
                  + d1rf6*w5 + Six*d1rw*f7*w5         &
                  + d1rf7*w6 + Seven*d1rw*f8*w6       &
                  + d1rf8*w7 + Eight*d1rw*f9*w7 + d1rf9*w8

!       The final value of term1 for 0 < omega < wcutoff is:

        term1 = t1 + t2t9 + t10

        d1sterm1 = d1st1 + d1st2t9 + d1st10
        d1rterm1 = d1rt1 + d1rt2t9 + d1rt10

!       The final value for the enhancement factor and its
!       derivatives is:

        Fx_wpbe = X * (term1 + term2 + term3 + term4 + term5)

        d1sfx = X * (d1sterm1 + d1sterm2 + d1sterm3    &
                              + d1sterm4 + d1sterm5)

        d1rfx = X * (d1rterm1 + d1rterm3 + d1rterm4 + d1rterm5)

      endif

      END SUBROUTINE wpbe_analy_erfc_approx_grad
