!
! Copyright (C) 2001-2009 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 pz_polarized (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and A. Zunger, PRB 23, 5048 (1981)
  !     spin-polarized energy and potential
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, b, c, d, gc, b1, b2
  parameter (a = 0.01555d0, b = - 0.0269d0, c = 0.0007d0, d = &
       - 0.0048d0, gc = - 0.0843d0, b1 = 1.3981d0, b2 = 0.2611d0)
  real(DP) :: lnrs, rs12, ox, dox
  REAL(DP), PARAMETER :: xcprefact = 0.022575584d0, pi34 = 0.6203504908994d0 
  ! REAL(DP) :: betha, etha, csi, prefact
  !
  if (rs.lt.1.0d0) then
     ! high density formula
     lnrs = log (rs)
     ec = a * lnrs + b + c * rs * lnrs + d * rs
     vc = a * lnrs + (b - a / 3.d0) + 2.d0 / 3.d0 * c * rs * lnrs + &
          (2.d0 * d-c) / 3.d0 * rs
  else
     ! interpolation formula
     rs12 = sqrt (rs)
     ox = 1.d0 + b1 * rs12 + b2 * rs
     dox = 1.d0 + 7.d0 / 6.d0 * b1 * rs12 + 4.d0 / 3.d0 * b2 * rs
     ec = gc / ox
     vc = ec * dox / ox
  endif
  !
!  IF ( lxc_rel ) THEN
!     betha = prefact * pi34 / rs
!     etha = DSQRT( 1 + betha**2 )
!     csi = betha + etha
!     prefact = 1.0D0 - (3.0D0/2.0D0) * ( (betha*etha - log(csi))/betha**2 )**2
!     ec = ec * prefact
!     vc = vc * prefact
!  ENDIF
  return
end subroutine pz_polarized
!
!-----------------------------------------------------------------------
subroutine pz_spin (rs, zeta, ec, vcup, vcdw)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, zeta, ec, vcup, vcdw
  !
  real(DP) :: ecu, vcu, ecp, vcp, fz, dfz
  real(DP) :: p43, third
  parameter (p43 = 4.0d0 / 3.d0, third = 1.d0 / 3.d0)
  !
  ! unpolarized part (Perdew-Zunger formula)
  call pz (rs, 1, ecu, vcu)
  ! polarization contribution
  call pz_polarized (rs, ecp, vcp)
  !
  fz = ( (1.0d0 + zeta) **p43 + (1.d0 - zeta) **p43 - 2.d0) / &
       (2.d0**p43 - 2.d0)
  dfz = p43 * ( (1.0d0 + zeta) **third- (1.d0 - zeta) **third) &
       / (2.d0**p43 - 2.d0)
  !
  ec = ecu + fz * (ecp - ecu)
  vcup = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * (1.d0 - zeta)
  vcdw = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * ( - 1.d0 - &
       zeta)
  !
  return
end subroutine pz_spin
!
!---------
SUBROUTINE vwn_spin(rs, zeta, ec, vcup, vcdw)

   USE kinds, ONLY: DP
   IMPLICIT NONE

   ! parameters:   e_c/para,    e_c/ferro,     alpha_c
   real(DP), parameter :: &
      A(3)  = (/ 0.0310907_dp, 0.01554535_dp, -0.01688686394039_dp /), &
      x0(3) = (/ -0.10498_dp, -0.32500_dp, -0.0047584_dp /), &
      b(3)  = (/3.72744_dp, 7.06042_dp, 1.13107_dp /), &
      c(3)  = (/ 12.9352_dp, 18.0578_dp, 13.0045_dp /),&
      Q(3)  = (/ 6.15199081975908_dp, 4.73092690956011_dp, 7.12310891781812_dp /), &
      tbQ(3) = (/ 1.21178334272806_dp, 2.98479352354082_dp, 0.31757762321188_dp /), &
      fx0(3) = (/ 12.5549141492_dp, 15.8687885_dp, 12.99914055888256_dp /), &
      bx0fx0(3) = (/ -0.03116760867894_dp, -0.14460061018521_dp, -0.00041403379428_dp /)
   ! N.B.: A is expressed in Hartree
   ! Q = sqrt(4*c - b^2)
   ! tbQ = 2*b/Q
   ! fx0 = X(x_0) = x_0^2 + b*x_0 + c
   ! bx0fx0 = b*x_0/X(x_0)

   real(DP), intent(in) :: rs, zeta
   real(DP), intent(out):: ec, vcup, vcdw

   ! local
   real(DP) :: zeta3, zeta4, trup, trdw, trup13, trdw13, fz, dfz, fzz4 
   real(DP) :: sqrtrs, ecP, ecF, ac, De, vcP, vcF, dac, dec1, dec2
   real(DP) :: cfz, cfz1, cfz2, iddfz0

   ! coefficients for f(z), df/dz, ddf/ddz(0)
   cfz = 2.0_dp**(4.0_dp/3.0_dp) - 2.0_dp
   cfz1 = 1.0_dp / cfz
   cfz2 = 4.0_dp/3.0_dp * cfz1
   iddfz0 = 9.0_dp / 8.0_dp *cfz
   sqrtrs = sqrt(rs)
   zeta3 = zeta**3
   zeta4 = zeta3*zeta
   trup = 1.0_dp + zeta
   trdw = 1.0_dp - zeta
   trup13 = trup**(1.0_dp/3.0_dp)
   trdw13 = trdw**(1.0_dp/3.0_dp)
   fz = cfz1 * (trup13*trup + trdw13*trdw - 2.0_dp)         ! f(zeta)
   dfz = cfz2 * (trup13 - trdw13)     ! d f / d zeta

   call padefit(sqrtrs, 1, ecP, vcP)            ! ecF = e_c Paramagnetic
   call padefit(sqrtrs, 2, ecF, vcF)            ! ecP = e_c Ferromagnetic
   call padefit(sqrtrs, 3, ac, dac)             ! ac = "spin stiffness"

   ac = ac * iddfz0
   dac = dac * iddfz0
   De = ecF - ecP - ac ! e_c[F] - e_c[P] - alpha_c/(ddf/ddz(z=0))
   fzz4 = fz * zeta4
   ec = ecP + ac * fz  + De * fzz4

   dec1 = vcP + dac*fz + (vcF - vcP - dac) * fzz4       ! e_c - (r_s/3)*(de_c/dr_s)
   dec2 = ac*dfz + De*(4.0_dp*zeta3*fz + zeta4*dfz)       ! de_c/dzeta

   ! v_c[s] = e_c - (r_s/3)*(de_c/dr_s) + [sign(s)-zeta]*(de_c/dzeta)
   vcup = dec1 + (1.0_dp - zeta)*dec2
   vcdw = dec1 - (1.0_dp + zeta)*dec2


   contains
   !---
   subroutine padefit(x, i, fit, dfit)
      !----
      ! implements formula [4.4] in:
      ! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)

      USE kinds, ONLY: DP
      implicit none

      ! input
      real(DP) :: x     ! x is sqrt(r_s) 
      integer :: i      ! i is the index of the fit

      ! output
      real(DP) :: fit, dfit
      ! Pade fit calculated in x and its derivative w.r.t. rho
      ! rs = inv((rho*)^(1/3)) = x^2
      ! fit  [eq. 4.4]
      ! dfit/drho = fit - (rs/3)*dfit/drs = ec - (x/6)*dfit/dx

      ! local
      real(DP) :: sqx, xx0, Qtxb, atg, fx
      real(DP) :: txb, txbfx, itxbQ

      sqx = x * x                          ! x^2 = r_s
      xx0 = x - x0(i)                      ! x - x_0
      Qtxb = Q(i) / (2.0_dp*x + b(i))      ! Q / (2x+b)
      atg = atan(Qtxb)                     ! tan^-1(Q/(2x+b))
      fx = sqx + b(i)*x + c(i)             ! X(x) = x^2 + b*x + c

      fit = A(i) * (  log(sqx/fx) + tbQ(i)*atg - &
            bx0fx0(i) * ( log(xx0*xx0/fx) + (tbQ(i) + 4.0_dp*x0(i)/Q(i)) * atg )  )

      txb = 2.0_dp*x + b(i)
      txbfx = txb / fx
      itxbQ = 1.0_dp / (txb*txb + Q(i)*Q(i))

      dfit = fit - A(i) / 3.0_dp + A(i)*x/6.0_dp * (  txbfx + 4.0_dp*b(i)*itxbQ + &
              bx0fx0(i) * ( 2.0_dp/xx0 - txbfx - 4.0_dp*(b(i)+2.0_dp*x0(i))*itxbQ )  )

   end subroutine

end subroutine

!-----------------------------------------------------------------------
subroutine pw_spin (rs, zeta, ec, vcup, vcdw)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, zeta, ec, vcup, vcdw
  ! xc parameters, unpolarised
  real(DP) :: a, a1, b1, b2, b3, b4, c0, c1, c2, c3, d0, d1
  parameter (a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = &
       3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, &
       c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0)
  ! xc parameters, polarised
  real(DP) :: ap, a1p, b1p, b2p, b3p, b4p, c0p, c1p, c2p, c3p, d0p, &
       d1p
  parameter (ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p &
       = 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p = &
       0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p &
       = 1.7697d0)
  ! xc parameters, antiferro
  real(DP) :: aa, a1a, b1a, b2a, b3a, b4a, c0a, c1a, c2a, c3a, d0a, &
       d1a
  parameter (aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = &
       3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a = &
       0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a &
       = 0.3969d0)
  real(DP) :: fz0
  parameter (fz0 = 1.709921d0)
  real(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz
  real(DP) :: om, dom, olog, epwc, vpwc
  real(DP) :: omp, domp, ologp, epwcp, vpwcp
  real(DP) :: oma, doma, ologa, alpha, vpwca
  !
  !     if(rs.lt.0.5d0) then
  ! high density formula (not implemented)
  !
  !     else if(rs.gt.100.d0) then
  ! low density formula  (not implemented)
  !
  !     else
  ! interpolation formula
  zeta2 = zeta * zeta
  zeta3 = zeta2 * zeta
  zeta4 = zeta3 * zeta
  rs12 = sqrt (rs)
  rs32 = rs * rs12
  rs2 = rs**2
  ! unpolarised
  om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2)
  dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 * rs32 &
       + 2.d0 * b4 * rs2)
  olog = log (1.d0 + 1.0d0 / om)
  epwc = - 2.d0 * a * (1.d0 + a1 * rs) * olog
  vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs) * olog - 2.d0 / &
       3.d0 * a * (1.d0 + a1 * rs) * dom / (om * (om + 1.d0) )
  ! polarized
  omp = 2.d0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2)
  domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs + 1.5d0 * b3p * &
       rs32 + 2.d0 * b4p * rs2)
  ologp = log (1.d0 + 1.0d0 / omp)
  epwcp = - 2.d0 * ap * (1.d0 + a1p * rs) * ologp
  vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs) * ologp - &
       2.d0 / 3.d0 * ap * (1.d0 + a1p * rs) * domp / (omp * (omp + 1.d0) &
       )
  ! antiferro
  oma = 2.d0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2)
  doma = 2.d0 * aa * (0.5d0 * b1a * rs12 + b2a * rs + 1.5d0 * b3a * &
       rs32 + 2.d0 * b4a * rs2)
  ologa = log (1.d0 + 1.0d0 / oma)
  alpha = 2.d0 * aa * (1.d0 + a1a * rs) * ologa
  vpwca = + 2.d0 * aa * (1.d0 + 2.d0 / 3.d0 * a1a * rs) * ologa + &
       2.d0 / 3.d0 * aa * (1.d0 + a1a * rs) * doma / (oma * (oma + 1.d0) &
       )
  !
  fz = ( (1.d0 + zeta) ** (4.d0 / 3.d0) + (1.d0 - zeta) ** (4.d0 / &
       3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0)
  dfz = ( (1.d0 + zeta) ** (1.d0 / 3.d0) - (1.d0 - zeta) ** (1.d0 / &
       3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) )
  !
  ec = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) &
       * fz * zeta4
  !
  vcup = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
       * fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
       zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
       * (1.d0 - zeta)

  vcdw = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
       * fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
       zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
       * (1.d0 + zeta)
  !      endif
  !
  return
end subroutine pw_spin
!
!-----------------------------------------------------------------------
subroutine pw_spin_vec (rs, zeta, evc, length)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  integer :: length
  real(DP) :: rs(length), zeta(length), evc(length,3)
  ! xc parameters, unpolarised
  real(DP) :: a, a1, b1, b2, b3, b4, c0, c1, c2, c3, d0, d1
  parameter (a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = &
       3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, &
       c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0)
  ! xc parameters, polarised
  real(DP) :: ap, a1p, b1p, b2p, b3p, b4p, c0p, c1p, c2p, c3p, d0p, &
       d1p
  parameter (ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p &
       = 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p = &
       0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p &
       = 1.7697d0)
  ! xc parameters, antiferro
  real(DP) :: aa, a1a, b1a, b2a, b3a, b4a, c0a, c1a, c2a, c3a, d0a, &
       d1a
  parameter (aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = &
       3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a = &
       0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a &
       = 0.3969d0)
  real(DP) :: fz0
  parameter (fz0 = 1.709921d0)
  real(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz
  real(DP) :: om, dom, olog, epwc, vpwc
  real(DP) :: omp, domp, ologp, epwcp, vpwcp
  real(DP) :: oma, doma, ologa, alpha, vpwca
  integer :: i
  !
  !     if(rs.lt.0.5d0) then
  ! high density formula (not implemented)
  !
  !     else if(rs.gt.100.d0) then
  ! low density formula  (not implemented)
  !
  !     else
  ! interpolation formula
  do i=1,length
     zeta2 = zeta(i) * zeta(i)
     zeta3 = zeta2 * zeta(i)
     zeta4 = zeta3 * zeta(i)
     rs12 = sqrt (rs(i))
     rs32 = rs(i) * rs12
     rs2 = rs(i)**2
     ! unpolarised
     om = 2.d0 * a * (b1 * rs12 + b2 * rs(i) + b3 * rs32 + b4 * rs2)
     dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs(i) + 1.5d0 * b3 * rs32 &
          + 2.d0 * b4 * rs2)
     olog = log (1.d0 + 1.0d0 / om)
     epwc = - 2.d0 * a * (1.d0 + a1 * rs(i)) * olog
     vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs(i)) * olog - 2.d0 / &
          3.d0 * a * (1.d0 + a1 * rs(i)) * dom / (om * (om + 1.d0) )
     ! polarized
     omp = 2.d0 * ap * (b1p * rs12 + b2p * rs(i) + b3p * rs32 + b4p * rs2)
     domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs(i) + 1.5d0 * b3p * &
          rs32 + 2.d0 * b4p * rs2)
     ologp = log (1.d0 + 1.0d0 / omp)
     epwcp = - 2.d0 * ap * (1.d0 + a1p * rs(i)) * ologp
     vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs(i)) * ologp - &
          2.d0 / 3.d0 * ap * (1.d0 + a1p * rs(i)) * domp / (omp * (omp + 1.d0) &
          )
     ! antiferro
     oma = 2.d0 * aa * (b1a * rs12 + b2a * rs(i) + b3a * rs32 + b4a * rs2)
     doma = 2.d0 * aa * (0.5d0 * b1a * rs12 + b2a * rs(i) + 1.5d0 * b3a * &
          rs32 + 2.d0 * b4a * rs2)
     ologa = log (1.d0 + 1.0d0 / oma)
     alpha = 2.d0 * aa * (1.d0 + a1a * rs(i)) * ologa
     vpwca = + 2.d0 * aa * (1.d0 + 2.d0 / 3.d0 * a1a * rs(i)) * ologa + &
          2.d0 / 3.d0 * aa * (1.d0 + a1a * rs(i)) * doma / (oma * (oma + 1.d0) &
          )
     !
     fz = ( (1.d0 + zeta(i)) ** (4.d0 / 3.d0) + (1.d0 - zeta(i)) ** (4.d0 / &
          3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0)
     dfz = ( (1.d0 + zeta(i)) ** (1.d0 / 3.d0) - (1.d0 - zeta(i)) ** (1.d0 / &
          3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) )
     !
     evc(i,3) = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) &
          * fz * zeta4
     !
     evc(i,1) = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
          * fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
          zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
          * (1.d0 - zeta(i))

     evc(i,2) = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
          * fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
          zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
          * (1.d0 + zeta(i))
  end do
  !      endif
  !
end subroutine pw_spin_vec
!
!-----------------------------------------------------------------------
subroutine becke88_spin (rho, grho, sx, v1x, v2x)
  !-----------------------------------------------------------------------
  ! Becke exchange: A.D. Becke, PRA 38, 3098 (1988) - Spin polarized case
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, grho, sx, v1x, v2x
  ! input: charge
  ! input: gradient
  ! output: the up and down energies
  ! output: first part of the potential
  ! output: the second part of the potential
  !
  real(DP) :: beta, third
  parameter (beta = 0.0042d0, third = 1.d0 / 3.d0)
  real(DP) :: rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee
  !
  rho13 = rho**third
  rho43 = rho13**4
  xs = 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 = grho / rho43 * ( - beta / dd)
  v1x = - (4.d0 / 3.d0) * xs2 * beta * rho13 * ee / dd2
  v2x = beta * (ee-dd) / (rho43 * dd2)
  !
  return
end subroutine becke88_spin
!
!-----------------------------------------------------------------------
subroutine perdew86_spin (rho, zeta, grho, sc, v1cup, v1cdw, v2c)
  !-----------------------------------------------------------------------
  ! Perdew gradient correction on correlation: PRB 33, 8822 (1986)
  ! spin-polarized case
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, 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, dd, ddd
  !
  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
  dd = (2.d0) **third * sqrt ( ( (1.d0 + zeta) * 0.5d0) ** (5.d0 / &
       3.d0) + ( (1.d0 - zeta) * 0.5d0) ** (5.d0 / 3.d0) )
  ddd = (2.d0) ** ( - 4.d0 / 3.d0) * 5.d0 * ( ( (1.d0 + zeta) &
       * 0.5d0) ** (2.d0 / 3.d0) - ( (1.d0 - zeta) * 0.5d0) ** (2.d0 / &
       3.d0) ) / (3.d0 * dd)
  ephi = exp ( - phi)
  sc = grho / rho43 * cn * ephi / dd
  v1cup = sc * ( (1.d0 + phi) * dcn / cn - ( (4.d0 / 3.d0) - &
       (7.d0 / 6.d0) * phi) / rho) - sc * ddd / dd * (1.d0 - zeta) &
       / rho
  v1cdw = sc * ( (1.d0 + phi) * dcn / cn - ( (4.d0 / 3.d0) - &
       (7.d0 / 6.d0) * phi) / rho) + sc * ddd / dd * (1.d0 + zeta) &
       / rho
  v2c = cn * ephi / rho43 * (2.d0 - phi) / dd
  !
  return
end subroutine perdew86_spin
!
!-----------------------------------------------------------------------
subroutine ggac_spin (rho, zeta, grho, sc, v1cup, v1cdw, v2c)
  !-----------------------------------------------------------------------
  ! Perdew-Wang GGA (PW91) correlation part - spin-polarized
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, 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, vcup, vcdw, t, expe, af, y, &
       xy, qy, s1, h0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, dh1, &
       ddh1, fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, dh0zup, &
       dh0zdw, dh1zup, dh1zdw
  !
  rs = pi34 / rho**third
  rs2 = rs * rs
  rs3 = rs * rs2
  call pw_spin (rs, zeta, ec, vcup, vcdw)
  kf = xkf / rs
  ks = xks * sqrt (kf)
  fz = 0.5d0 * ( (1.d0 + zeta) ** (2.d0 / 3.d0) + (1.d0 - zeta) ** ( &
       2.d0 / 3.d0) )
  fz2 = fz * fz
  fz3 = fz2 * fz
  fz4 = fz3 * fz
  dfz = ( (1.d0 + zeta) ** ( - 1.d0 / 3.d0) - (1.d0 - zeta) ** ( - &
       1.d0 / 3.d0) ) / 3.d0
  t = sqrt (grho) / (2.d0 * fz * ks * rho)
  expe = exp ( - 2.d0 * al * ec / (fz3 * be * be) )
  af = 2.d0 * al / be * (1.d0 / (expe-1.d0) )
  bfup = expe * (vcup - ec) / fz3
  bfdw = expe * (vcdw - ec) / fz3
  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 = fz3 * be * be / (2.d0 * al) * log (s1)
  dh0up = be * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfup / be-7.d0 / 3.d0) )
  dh0dw = be * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfdw / be-7.d0 / 3.d0) )
  dh0zup = (3.d0 * h0 / fz - be * t * t * fz2 / s1 * (2.d0 * xy - &
       qy * (3.d0 * af * expe * ec / fz3 / be+2.d0) ) ) * dfz * (1.d0 - &
       zeta)
  dh0zdw = - (3.d0 * h0 / fz - be * t * t * fz3 / s1 * (2.d0 * xy - &
       qy * (3.d0 * af * expe * ec / fz3 / be+2.d0) ) ) * dfz * (1.d0 + &
       zeta)
  ddh0 = be * fz / (2.d0 * ks * ks * rho) * (xy - qy) / s1
  ee = - 100.d0 * fz4 * (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) * fz3 * t * t * exp (ee)
  dh1 = - third * (h1 * (7.d0 + 8.d0 * ee) + fz3 * nu * t * t * exp &
       (ee) * dcn)
  ddh1 = 2.d0 * h1 * (1.d0 + ee) * rho / grho
  dh1zup = (1.d0 - zeta) * dfz * h1 * (1.d0 + 2.d0 * ee / fz)
  dh1zdw = - (1.d0 + zeta) * dfz * h1 * (1.d0 + 2.d0 * ee / fz)
  sc = rho * (h0 + h1)
  v1cup = h0 + h1 + dh0up + dh1 + dh0zup + dh1zup
  v1cdw = h0 + h1 + dh0up + dh1 + dh0zdw + dh1zdw
  v2c = ddh0 + ddh1
  return
end subroutine ggac_spin
!
!---------------------------------------------------------------
subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
  !---------------------------------------------------------------
  !
  ! PBE correlation (without LDA part) - spin-polarized
  ! iflag = 1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
  ! iflag = 2: J.P.Perdew et al., PRL 100, 136406 (2008)
  !
  USE kinds, ONLY : DP
  implicit none
  integer, intent(in) :: iflag
  real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c
  real(DP) :: ga, be(2)
  parameter (ga = 0.031091d0)
  data be / 0.066725d0 ,  0.046d0 /
  real(DP) :: third, pi34, xkf, xks
  parameter (third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0)
  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, vcup, vcdw, t, expe, af, y, xy, qy, &
       s1, h0, ddh0
  real(DP) :: fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, &
       dh0zup, dh0zdw
  !
  rs = pi34 / rho**third
  call pw_spin (rs, zeta, ec, vcup, vcdw)
  kf = xkf / rs
  ks = xks * sqrt (kf)
  fz = 0.5d0 * ( (1.d0 + zeta) ** (2.d0 / 3.d0) + (1.d0 - zeta) ** ( &
       2.d0 / 3.d0) )
  fz2 = fz * fz
  fz3 = fz2 * fz
  fz4 = fz3 * fz
  dfz = ( (1.d0 + zeta) ** ( - 1.d0 / 3.d0) - (1.d0 - zeta) ** ( - &
       1.d0 / 3.d0) ) / 3.d0
  t = sqrt (grho) / (2.d0 * fz * ks * rho)
  expe = exp ( - ec / (fz3 * ga) )
  af = be(iflag) / ga * (1.d0 / (expe-1.d0) )
  bfup = expe * (vcup - ec) / fz3
  bfdw = expe * (vcdw - ec) / fz3
  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 = fz3 * ga * log (s1)
  dh0up = be(iflag) * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfup / be(iflag)-7.d0 / 3.d0) )
  dh0dw = be(iflag) * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfdw / be(iflag)-7.d0 / 3.d0) )
  dh0zup = (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
  qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 - zeta)
  dh0zdw = - (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
  qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 + zeta)

  ddh0 = be(iflag) * fz / (2.d0 * ks * ks * rho) * (xy - qy) / s1
  sc = rho * h0
  v1cup = h0 + dh0up + dh0zup
  v1cdw = h0 + dh0dw + dh0zdw
  v2c = ddh0
  return
end subroutine pbec_spin
!
!-----------------------------------------------------------------------
subroutine slater_spin (rho, zeta, ex, vxup, vxdw)
  !-----------------------------------------------------------------------
  !     Slater exchange with alpha=2/3, spin-polarized case
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, zeta, ex, vxup, vxdw
  real(DP) :: f, alpha, third, p43
  parameter (f = - 1.10783814957303361d0, alpha = 2.0d0 / 3.0d0)
  ! f = -9/8*(3/pi)^(1/3)
  parameter (third = 1.d0 / 3.d0, p43 = 4.d0 / 3.d0)
  real(DP) :: exup, exdw, rho13
  !
  rho13 = ( (1.d0 + zeta) * rho) **third
  exup = f * alpha * rho13
  vxup = p43 * f * alpha * rho13
  rho13 = ( (1.d0 - zeta) * rho) **third
  exdw = f * alpha * rho13
  vxdw = p43 * f * alpha * rho13
  ex = 0.5d0 * ( (1.d0 + zeta) * exup + (1.d0 - zeta) * exdw)
  !
  return
end subroutine slater_spin

!-----------------------------------------------------------------------
subroutine slater_spin_vec(rho, zeta, evx, length)
  !-----------------------------------------------------------------------
  !     Slater exchange with alpha=2/3, spin-polarized case
  !
  USE kinds, ONLY : DP
  implicit none
  integer  :: length
  real(DP) :: rho(length), zeta(length), evx(length,3)
  real(DP) :: f, alpha, third, p43
  parameter (f = - 1.10783814957303361d0, alpha = 2.0d0 / 3.0d0)
  ! f = -9/8*(3/pi)^(1/3)
  parameter (third = 1.d0 / 3.d0, p43 = 4.d0 / 3.d0)
  real(DP) :: exup(length), exdw(length), rho13(length)
  !
  rho13 = ( (1.d0 + zeta) * rho) **third
  exup = f * alpha * rho13
  evx(:,1) = p43 * f * alpha * rho13
  rho13 = ( (1.d0 - zeta) * rho) **third
  exdw = f * alpha * rho13
  evx(:,2) = p43 * f * alpha * rho13
  evx(:,3) = 0.5d0 * ( (1.d0 + zeta) * exup + (1.d0 - zeta) * exdw)
  !
end subroutine slater_spin_vec

!-----------------------------------------------------------------------
SUBROUTINE slater_rxc_spin ( rho, Z, ex, vxup, vxdw )
  !-----------------------------------------------------------------------
  !     Slater exchange with alpha=2/3, relativistic exchange case
  !
  USE kinds, ONLY : DP
  USE constants, ONLY : pi
  IMPLICIT none
  real (DP):: rho, ex, vxup, vxdw
  !
  real(DP), PARAMETER :: ZERO=0.D0, ONE=1.D0, PFIVE=.5D0, &
       OPF=1.5D0, C014=0.014D0
  real (DP):: rs, trd, ftrd, tftm, a0, alp, z, fz, fzp, vxp, xp, &
       beta, sb, alb, vxf, exf

  TRD = ONE/3.d0
  FTRD = 4.d0*TRD
  TFTM = 2**FTRD-2.d0
  A0 = (4.d0/(9.d0*PI))**TRD
  
  !      X-alpha parameter:
  ALP = 2.d0 * TRD
  
  IF (rho <=  ZERO) THEN
     EX = ZERO
     vxup  = ZERO
     vxdw  = ZERO
     RETURN
  ELSE
     FZ = ((1.d0+Z)**FTRD+(1.d0-Z)**FTRD-2.d0)/TFTM
     FZP = FTRD*((1.d0+Z)**TRD-(1.d0-Z)**TRD)/TFTM
  ENDIF
  RS = (3.d0 / (4.d0*PI*rho) )**TRD
  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.d0**TRD*VXP
  EXF = 2.d0**TRD*XP
  vxup  = VXP + FZ*(VXF-VXP) + (1.d0-Z)*FZP*(EXF-XP)
  vxdw  = VXP + FZ*(VXF-VXP) - (1.d0+Z)*FZP*(EXF-XP)
  EX    = XP + FZ*(EXF-XP)
       
END SUBROUTINE slater_rxc_spin


!-----------------------------------------------------------------------
subroutine slater1_spin (rho, zeta, ex, vxup, vxdw)
  !-----------------------------------------------------------------------
  !     Slater exchange with alpha=2/3, spin-polarized case
  !
  use kinds, only: dp
  implicit none
  real(DP) :: rho, zeta, ex, vxup, vxdw
  real(DP), parameter :: f = - 1.10783814957303361d0, alpha = 1.0d0, &
       third = 1.d0 / 3.d0, p43 = 4.d0 / 3.d0
  ! f = -9/8*(3/pi)^(1/3)
  real(DP) :: exup, exdw, rho13
  !
  rho13 = ( (1.d0 + zeta) * rho) **third
  exup = f * alpha * rho13
  vxup = p43 * f * alpha * rho13
  rho13 = ( (1.d0 - zeta) * rho) **third
  exdw = f * alpha * rho13
  vxdw = p43 * f * alpha * rho13
  ex = 0.5d0 * ( (1.d0 + zeta) * exup + (1.d0 - zeta) * exdw)
  !
  return
end subroutine slater1_spin
!
!-----------------------------------------------------------------------
function dpz_polarized (rs, iflg)
  !-----------------------------------------------------------------------
  !  derivative of the correlation potential with respect to local density
  !  Perdew and Zunger parameterization of the Ceperley-Alder functional
  !  spin-polarized case
  !
  USE kinds, only : DP
  USE constants, ONLY : pi, fpi
  !
  implicit none
  !
  real(DP), intent (in) :: rs
  integer, intent(in) :: iflg
  real(DP) :: dpz_polarized
  !
  !  local variables
  !  a,b,c,d,gc,b1,b2 are the parameters defining the functional
  !
  real(DP), parameter :: a = 0.01555d0, b = -0.0269d0, c = 0.0007d0, &
       d = -0.0048d0, gc = -0.0843d0, b1 = 1.3981d0, b2 = 0.2611d0,&
       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_polarized = - fpi * rs**4.d0 / 9.d0 * dmrs
  return
  !
end function dpz_polarized